PART 1: THE MANUSCRIPT

Short running title

Telomere length attrition predicts health span

Authors & Affiliations:

Luise A. Seeker1@, Sarah L. Underwood2, Rachael V. Wilbourn2, Jennifer Fairlie2, Hannah Froy2@@, Rebecca Holland2, Joanna J. Ilska1, 3, Androniki Psifidi3,4, Ainsley Bagnall5, Bruce Whitelaw3, Mike Coffey1, Georgios Banos1, 3 & Daniel H. Nussey2

1 Animal & Veterinary Sciences, SRUC, Roslin Institute Building, Easter Bush, Midlothian, UK 2 Institute of Evolutionary Biology, School of Biological Sciences, University of Edinburgh, UK 3 The Roslin Institute and Royal (Dick) School of Veterinary Studies, University of Edinburgh, Easter Bush, Midlothian, UK 4 Royal Veterinary College, University of London, Hatfield, UK 5 SRUC Crichton Royal Farm, Glencaple Road, Dumfries, UK

Corresponding author:

Luise A. Seeker MRC Centre for Regenerative Medicine, The University of Edinburgh, Edinburgh BioQuarter, 5 Little France Drive, Edinburgh EH16 4UU Luise.seeker[AT]ed.ac.uk 07591133397 Present author addresses: @ MRC Centre for Regenerative Medicine, University of Edinburgh, Edinburgh, UK @@ Centre for Biodiversity Dynamics, NTNU Norwegian University of Science and Technology, Trondheim, Norway

Abstract

Telomere length measured in blood cells is predictive of subsequent adult health and survival across a range of vertebrate species. However, we currently do not know whether such associations result from among-individual differences in telomere length determined genetically or by environmental factors early in life, or from differences in the rate of telomere attrition over the course of life. Here, we measured relative leukocyte telomere length (RLTL) multiple times across the entire lifespan of dairy cattle in a research population that is closely monitored for health and milk production and where individuals are only culled in response to health issues and less due to poor milk production than on purely commercial farms. Our results clearly show that the average amount of telomere attrition over an individual’s life, not their average or early life telomere length predicted when an individual was culled. Within-individual telomere length attrition could reflect environmental or physiological insults which may accumulate to predict individual health-span. We also show that animals with more telomere attrition in their first year of life were culled at a younger age, indicating that early life stressors may have a prolonged effect on adult life.

Keywords:

Telomeres, telomere length, longitudinal, lifespan, health-span, qPCR

Introduction

Telomeres are repetitive DNA sequences that cap the ends of eukaryote linear chromosomes (Blackburn & Gall, 1978; De Lange, 2005). They shorten with the number of cell divisions in vitro as well as in response to oxidative stress and critically short telomeres trigger a DNA damage response that leads to replicative senescence or apoptosis (Harley, Futcher, & Greider, 1990; Olovnikov, 1973; Watson, 1972). In the last decade or so, measures of average telomere length (TL) taken from blood samples have emerged as an exciting biomarker of health across disciplines including biomedicine, epidemiology, ecology and evolutionary biology (Aviv & Shay, 2018; Harrington & Pucci, 2018; Wilbourn et al., 2018). Considerable among- and within-individual variation in TL has been observed, with a general pattern of rapid telomere attrition during early life and a plateau or slower decline thereafter (Aubert & Lansdorp, 2008; Baerlocher, Rice, Vulto, & Lansdorp, 2007). Both genetic and environmental factors, particularly those associated with physiological stress, predict TL in humans and other vertebrates (Angelier, Costantini, Blévin, & Chastel, 2017; Muhammad Asghar et al., 2015; Cherkas et al., 2006; Dugdale, Richardson, & Richardson, 2018; Epel et al., 2004). TL has also been repeatedly associated with health outcomes and subsequent survival in a variety of species, particularly humans and birds (Boonekamp, Simons, Hemerik, & Verhulst, 2013; Wilbourn et al., 2018) and also experimentally elongated TL in mice were associated with a survival advantage (Muñoz-lorente, Cano-martin, & Blasco, 2019). However, a major outstanding question remains to what degree associations between TL and health arise from constitutive differences in TL among individuals set by genes or early life conditions, or from the pattern of within-individual change in TL across individuals’ lives. Estimates of the individual consistency of TL over time in both human and avian literature vary considerably among studies. Some studies report very high correlations within individuals across follow-up measurements (r > 0.8; (Benetos et al., 2013; Boonekamp, Mulder, Salomons, Dijkstra, & Verhulst, 2014)) and high heritability of TL (>0.7; (Atema et al., 2015; Broer et al., 2013)), and have shown in humans that the rank order in TL among individuals remains relatively consistent during adult life (Benetos et al., 2013). This implies most of the variation in blood cell TL occurs at the among-individual level and is predominantly determined by genetics and early-life environment (Atema et al., 2015; Benetos et al., 2013). In stark contrast to this, a growing body of literature studying both humans and non-human vertebrates reports much lower individual repeatability and heritability (Dugdale & Richardson, 2018; Voillemot et al., 2012). These studies demonstrate very high levels of within-individual variation in TL and show that changes in TL across consecutive measurements are highly dynamic (Fairlie et al., 2015; Farzaneh-Far et al., 2010). Although on average TL attrition over time tends to be the norm, a growing number of studies have found that a substantial proportion of individuals show telomere elongation over time and it is now clear this cannot be ascribed solely to measurement error (Bateson & Nettle, 2016; Spurgin et al., 2018). A recent cross-sectional study has shown that amongst different avian and mammal species, those with a fast telomere attrition rate have a shorter life span but without investigating if telomere attrition rate within individuals is predictive of their life span within the population as well (Whittemore, Vera, Martínez-nevado, Sanpera, & Blasco, 2019). To date, few studies have directly compared the power of an individual’s average or early life TL versus the pattern of within-individual change in TL to explain variation in measured associated with health and fitness. In Seychelle warblers and Alpine swifts both telomere length and telomere attrition rate were associated with survival (Barrett, Burke, Hammers, Komdeur, & Richardson, 2013; Bize, Criscuolo, Metcalfe, Nasir, & Monaghan, 2009) while in jackdaw nestlings telomere length was not predictive of lifespan but early life telomere attrition- accelerated by early life stress- was (Boonekamp, Mulder, Salomons, Dijkstra, & Verhulst, 2014). Contrary, in a study on elderly humans, there was an association between telomere length and longevity, but not between telomere attrition and life span (Yuan, Kronström, Hellenius, Cederholm, & Xu, 2018). To our knowledge no studies so far has investigated the relationship between telomere attrition and life span with data spanning the entire typical lifespan of the studied individuals. Here, we use an exceptionally detailed longitudinal study of dairy cattle to show that a life-long tendency for greater telomere attrition, rather than average TL, predicts productive lifespan, a metric closely linked to long-term health.

We used samples and data collected as part of the long-term study of Holstein Friesian dairy cattle kept at the Crichton Royal Research Farm in Dumfries, Scotland (Veerkamp, Simm, & Oldham, 1994). This herd, consisting of around 200 milking cows plus their calves and replacement heifers, has been regularly monitored since 1973 for a broad range of measurements, such as body weight, feed intake, signs of disease (health events), milk yield, productive lifespan and reasons for culling (Veerkamp et al., 1994). The cows are separated into two distinct genetic groups (a group that has been strongly selected for high milk fat and protein yield, and a control group) and are randomly allocated to two different diets (a high forage, low energy diet vs. a low forage, high energy diet group) (Veerkamp et al., 1994). Blood samples are routinely collected from members of the herd, initially within 15 days after birth and then approximately annually thereafter. Productive lifespan, which is the age of an individual at culling, is recorded for every animal together with a reason for culling. In this and other commercial livestock herds, culling is a response to diseases such as reproductive problems, lameness and mastitis. Although productive lifespan in livestock is clearly not a comparable measure to longevity in study of humans or wild vertebrates, we argue that because age at culling is tightly linked to health, it can be considered as a proxy for health-span, over and above its obvious importance as a metric of agricultural productivity. We have previously reported for this herd that TL is moderately repeatable and heritable, but that individual slopes relating TL to age vary considerably among individuals (Seeker et al., 2018; Seeker et al., 2018). We have also reported that TL at the age of one year predicts productive lifespan in this population of dairy cattle (Seeker et al., 2018). Here, we test the relationships between lifetime average TL as well as lifetime change in TL with productive lifespan.

Materials and Methods

Animal population and data collection

At the Crichton Royal Farm, 200 milking cows plus their calves and replacement heifers are kept at any time. One half of the milking cows belong to a genetic line that has been selected for high milk protein and fat yield (S), while the other half is deliberately maintained on a UK average productivity level (C). Calves and heifers of both genetic lines are kept together. After first calving all cows are randomly allocated to a high forage (HF) or low forage (LF) diet. The LF diet is energy richer than the HF diet and whilst the LF cows are housed continuously, the HF cows graze over the summer months. All cows are milked three times daily and milk yield is recorded. In the present study, these measurements were used to calculate an average milk production in kg per cow including all started lactations and it is referred to this below as “average lifetime milk production” (Figure S1). Every day cows leave the milking parlour over a pressure plate which detects signs of lameness. Behaviour and health events are documented after visual detection by farm workers (Figure S2). At the end of the animal’s life its productive lifespan and a reason for culling are recorded. Productive lifespan is the time from birth to culling in days and is a proxy for the health span of the animal, because all animals that remain healthy enough to generate profit for the farmer remain in the herd. The most frequent reasons for culling were reproductive problems, mastitis, lameness and injuries caused by accidents (Figure S3). Routine blood sampling takes place initially shortly after birth and then annually in spring. Because of this sampling routine and because calves are being born all year round, age at sampling and sampling intervals vary for adult animals (Figure S1). Further information about the animal population can be found in Supplementary File 2.

DNA extraction and qPCR

DNA from whole blood samples was extracted with the DNeasy Blood and Tissue spin column kit (QIAGEN) and telomere length was measured by qPCR as previously described (Seeker et al., 2016; Seeker et al., 2018; Seeker et al., 2018; Seeker, Holland, Psifidi, Banos, & Nussey, 2015). The repeatability of the assay (see Supplementary File 2 for how repeatability was calculated) was 80% and therefore delivers robust results for further interpretation (Nettle, Seeker, Nussey, Froy, & Bateson, 2019). A full description of our DNA extraction and qPCR protocols including quality control steps can also be found in Supplementary File 2.

Statistical analysis

We have shown before that our RLTL data are significantly affected by qPCR plate and qPCR row (Seeker et al., 2016; Seeker et al., 2018). To account for those known sources of measurement error, we used the residuals of a linear model that corrected all RLTL measurements for qPCR plate and row, by fitting plate and row as fixed factors in the model. These residual RLTL measures were used in all subsequent calculations and models of telomere dynamics. RLTL change was calculated as the difference between two subsequent adjusted RLTL measurements within individuals (RLTL change = RLTLt - RLTLt-1). To investigate which factors affect the direction and amount of RLTL change, a linear mixed model was fitted with RLTL change as response variable and animal identity as random effect. To capture the experimental design, feed group and genetic group were added as fixed effects to all models. Also, the time difference between consecutive samplings in days was fitted as a covariate, becasue sampling intervals varied within and between animals. The following additional factors were tested for their relationship with RLTL change: birth year of the animal (fixed effect), year of sampling (fixed effect), age at sampling (at time t), and the occurrence of a health event within two weeks before or after sampling (at time t).

Becasue of our study design, the first two sample years only included RLTL change in young animals which has previously been shown to follow a different pattern. To exclude that a year effect was driven by the difference in juvenile and adult RLTL dynamics, the full model was repeated while excluding the first two sample years.

Birth year and the occurrence of a heath event were not statistically significant (p>0.05) and therefore eliminated from the model. In the reduced model, age was first fitted as a covariate (age in years) than as a fixed effect (factor with two levels: younger vs. older than 1 year) and model fit was compared using AIC values. Because we wanted to investigate if milk productivity was associated with RLTL change, the best-fitting model of RLTL change was repeated for 918 change measurements of 253 animals for which milk productivity measurements were available and average lifetime milk production re-scaled by dividing it by 1000 was fitted as additional covariate.

To investigate the association between RLTL change and productive lifespan we first focussed on RLTL change within the first year of life by only considering two RLTL measurements per animal: The first was taken shortly after birth and the second at the approximate age of one year. The change between those measurements was noramlly distributed and tested as an explanatory variable in a Cox proportional hazard model of productive lifespan.

Our data shows similar to other longitudinal studies with more than two measurements no constant attrition, maintenance or elongation, but more complex dynamics with short term changes in both directions. Example individual RLTL dynamics are shown for all cows with at least seven RLTL measurements in our dataset in figure S4. The complexity of telomere dynamics (which cannot be explained by measurement error alone) means that it is impossible to compare lifelong telomere dynamics by simply comparing slopes. To still be able to test how life-long telomere dynamics are associated with productive lifespan, we calculated three average metrics of RLTL dynamics over every individual’s lifetime: Firstly, the average of all their RLTL measures (“mean RLTL”) which investigates if animals with on average longer RLTL have a survival advantage, secondly, the average of all their RLTL change measures (“mean RLTL change”) and thirdly, the average of all their absolute RLTL change measures (“mean absolute RLTL change”). While mean RLTL change captures the direction and magnitude of RLTL changes, mean absolute RLTL change describes just the magnitude of change without considering its direction because we were interested to investigate if more change in either direction may be correlated with adverse effects. Figure S5 visualises the reasoning behind calculating these three measures of lifetime telomere length dynamics which are not surprisingly moderately correlated with one another (r ranged from -0.53 to 0.30, Figure S6). We used Cox proportional hazard models of productive lifespan that also included mean milk production as covariate and included mean RLTL, mean RLTL change and mean absolute RLTL change as explanatory variables first separately and then together to test their association with productive lifespan first separately, then while accounting for the effect of the other two measures. We have reported before that on average RLTL shortens in this population of dairy cattle within the first year of life while remaining relatively stable thereafter at the population level (Seeker et al., 2018; Seeker, et al., 2018). We were concerned that the initial RLTL attrition may contribute proportionally more to overall RLTL shortening in short-lived animals with less follow up samples which may explain a correlation of more telomere attrition with a shorter productive lifespan. To ensure the relationship was not entirely driven by this correlation, we tested all models again while excluding all measurements that were taken shortly after birth and then also while excluding all animals with less than three RLTL measurements. Lastly, we tested if the previously reported effect of RLTL at the age of one year on productive lifespan (Seeker, Ilska, Psifidi, Wilbourn, Underwood, Fairlie, Holland, Froy, Salvo-Chirnside, et al., 2018) remained statistically significant when tested in a Cox proportional hazard model together with milk productivity and mean RLTL change. All statistical analyses were performed in R studio with R 3.1.3. (R Core Team, 2014). Mixed-effects models were implemented using the ‘lme4’ library (Bates, Mächler, Bolker, & Walker, 2015), while Cox proportional hazard models were implemented using the library survival (Therneau, 2015) and figures were generated with the library ‘ggplot2’ (Wickham, 2009).

Results

We used blood samples collected between 2008 and 2014 to measure relative leukocyte telomere length (RLTL) by monoplex qPCR in 1,325 samples from 305 individuals. On average 4.3 (range: 2-8) telomere measurements were made of each individual, including the first measurement within 15 days of birth and a variable number of subsequent measures (Figure S7 A). RLTL measures, adjusted for qPCR plate and row (see Methods), were approximately normally distributed (Figure S7 B).

From these longitudinal data we calculated 1,020 measures of change in adjusted RLTL across subsequent measurements of individual animals. RLTL change also showed approximately a normal distribution with a mean statistically significantly smaller than zero (Figure S7 C; t = -5.844, df = 1019, p-value <0.001) indicating that telomere shortening overall predominated.

Animals varied in the amount and direction of RLTL change across consecutive measurements, with a fairly even proportion of individuals increasing (43.5%) and decreasing (56.5%) in RLTL over time.

Modelling

Factors affecting telomere change

We next ran a series of models to test whether known individual, genetic and environmental variables could explain variation in RLTL change. Only sample year was significant in these models; age in years, genetic group, feed group, birth year, the time difference between samplings in days, and the occurrence of a health event within two weeks of sampling were not significant (Table S1). However, when age was fitted a factor with two levels for animals younger or oler than 1 year, it became clear that younger anumals show significantly more telomere length attrition (Tables S2 and S3). Fitting age as a covariate (age in years) and not a two-level factor lead to a slightly better although not statistically significant model fit and was therefore realised in all following models (AICs: -601.82 vs. -600.99).

We have previously shown that, on average across this population, RLTL declined over the first year of life but showed no systematic change with age thereafter (Seeker, Ilska, Psifidi, Wilbourn, Underwood, Fairlie, Holland, Froy, Salvo-Chirnside, et al., 2018; Seeker, Ilska, Psifidi, Wilbourn, Underwood, Fairlie, Holland, Froy, Bagnall, et al., 2018). Consistent with this, we found that average RLTL change across consecutive measurements was only significantly negative (indicating a tendency for attrition over time) when the first measurement was made close to birth and the follow up measurement at the age of around 1 year (Figure 1A, estimate= -0.115, SE= 0.01, t= -11.512, p<0.001; Table S3).

Average lifetime milk productivity was not statistically associated with RLTL change, when tested in a subset of animals (253 animals with 918 RLTL change measurements) that had milk productivity measurements available (Table S2).

Consecutive RLTL measurements made on the same individual were overall moderately positively and significantly correlated (r=0.38, p<0.001; Figure 1B), supporting our previously reported moderate and significant individual repeatability of RLTL (Seeker, Ilska, Psifidi, Wilbourn, Underwood, Fairlie, Holland, Froy, Bagnall, et al., 2018). We next tested whether measured RLTL change predicted the productive lifespan of cattle. Of all dead cows (N=244) the vast majority (N=241) had survived to their first lactation, but there was considerable variation in productive lifespan beyond this point (Figure S7 D). We used a Cox proportional hazard model to test whether early life RLTL change between two samples, one taken shortly after birth and the next at an approximate age of one year, predicted productive lifespan and found that greater early life RLTL attrition was associated with a shorter productive lifespan (N= 291, coefficient= -1.141, SE= 0.391, p= 0.004, Figure 2, Table S4). When each measure of lifetime RLTL dynamics was tested in a separate cox proportional hazard model of productive lifespan, mean RLTL did not significantly predict productive lifespan (coefficient= 0.341, SE= 0.591, p=0.564, Figure 2B , Table S4), but both mean RLTL change (coefficient= -5.209, SE= 0.845, p<0.001; Figure 2D, Table S4) and mean absolute RLTL change (coefficient= 2.939, SE= 0.970, p=0.002; Figure 2C, Table S4) were significantly associated with productive lifespan. When all three measures of lifetime RLTL dynamics were included in the same model only mean RLTL change remained significant (coefficient = -4.758, SE= 1.018, p<0.001; Table S4). This implies that the relationship between productive lifespan and mean absolute RLTL change was largely due to covariance with mean RLTL change. Thus, individuals that experienced greater telomere attrition over their lifetimes had a shorter productive lifespan and direction of RLTL change (rather than simply absolute magnitude) was an important aspect of this relationship. To visualise the relationship between RLTL change measurements and productive lifespan using Kaplan-Meier plots (Figure 2), continuous RLTL measures were transformed to a discrete scale by grouping them into tertiles (Figure S8). Cox proportional hazard models based on these tertile groupings of RLTL measures showed similar results to those reported above (Figure 2). The relationship between mean RLTL change and productive lifespan was robust to the inclusion of milk production (a physiological stressor, which is positively associated with productive lifespan as cows with a low milk yield are more likely to be culled at a younger age) in the model (Table S5). If RLTL declines most within the first year of life (as Figure 1A indicates), the association between change in RLTL and productive lifespan may be driven by the fact that this initial decline contributes relatively more to estimates of mean RLTL change in shorter-lived individuals than in animals that have more follow up samples with more moderate change measures available. We therefore repeated the analysis excluding the early life RLTL measurements and found that mean RLTL change still predicted productive lifespan (N=253, coefficient= -5.056, SE= 1.315, p<0.001, Table S6). The relationship between mean RLTL change and productive lifespan also remained statistically significant, when animals with less than 3 samples were excluded from the analysis (N= 213, coefficient= -3.47, SE= 1.34, Wald test = 6.7 on 1 df, p= 0.01) and when RLTL at 1 year of age (which we previously found to predict productive lifespan (Seeker, Ilska, Psifidi, Wilbourn, Underwood, Fairlie, Holland, Froy, Salvo-Chirnside, et al., 2018)) was included in the models (Table S7).

Mean rate of RLTL change and mean rate of absolute RLTL change correlate well with mean RLTL change and mean absolute RLTL change. THerefore, they do not seem to add much to the analysis. Only mean RLTL, mean RLTL change and mean absolute RLTL change are used.

Discussion

We have presented, to our knowledge, the first demonstration in any vertebrate that lifelong variation in telomere attrition rather than variation in constitutive individual differences in average telomere length predict health outcomes and lifespan. While there is mounting evidence that telomere length predicts mortality, health and life history in humans as well as birds and non-human mammals (Bakaysa et al., 2007; Bichet et al., 2019; Boonekamp et al., 2013; Cawthon, Smith, O’Brien, Sivatchenko, & Kerber, 2003; Fairlie et al., 2015; Mark Haussmann, Winkler, & Vleck, 2005; Heidinger et al., 2012; Nettle, Monaghan, Boner, Gillespie, & Bateson, 2013; Voillemot et al., 2012; Wilbourn et al., 2018), very few studies have been able to accumulate long-term longitudinal data capable of differentiating the role of among- and within-individual variation in TL to such relationships. Our data support the contention that within-individual directional change over time in TL is more important than stable among-individual differences in predicting overall health. Individuals in our study which showed a greater propensity to shorten RLTL across sampling points were culled earlier, mainly due to health issues described above. There was no relationship between productive lifespan and an individual’s average RLTL. Future studies will show how well our results generalise to other systems as telomere biology is variable amongst species. Cattle telomere biology seems to be similar to other ruminants, horses, zebras, tapirs some whales and primates including humans in that they have relatively short telomeres and a tight regulation of telomerase expression (Gomes et al., 2011). If our results generalise to some of those other systems and contexts, they have important implications for the utility of TL as a biomarker of health and fitness, lending support to the idea that TL change is an indirect marker reflecting past physiological insults and stress rather than an indicator of constitutive or genetically-based robustness to life’s challenges. Our data also highlight the importance of collecting longitudinal telomere measurements, by showing that in some species it is within-individual change over time in TL that carries the important biological signal.

Our study animals varied considerably in the magnitude and direction of RLTL change across consecutive sampling points. This is in accordance with other longitudinal studies that have reported a wide variation in TL change, alongside evidence that a large proportion of individuals show telomere lengthening over time (Bize et al., 2009; Chen et al., 2011; Farzaneh-Far et al., 2010; Gardner et al., 2005; Huzen et al., 2014; Kark, Goldberger, Kimura, Sinnreich, & Aviv, 2012; Nordfjäll et al., 2009; Shalev et al., 2012; Steenstrup et al., 2013; Svenson et al., 2011). Whilst the vast majority of studies linking telomere length with health and fitness remain cross-sectional, a growing number of longitudinal studies have found that individuals that lose more TL across samples have reduced survival prospects (Barrett et al., 2013; Bize et al., 2009; Boonekamp et al., 2014; Goglin et al., 2016). However, these studies tend to involve two samples taken per individual or a follow-up period that represents only a small fraction of the total lifespan of the study species. A lifelong study of TL and lifespan in zebra finches found that TL measured in early life was a better predictor of subsequent survival than measures from later life, but the role of change in telomere length across sampling points was not investigated (Heidinger et al., 2012). Here we show that, despite TL being moderately consistent across the lifetimes of individuals, there remains considerable within-individual variation and the pattern of change in TL over an individual’s life is highly dynamic. Short-term environmental fluctuations impacting TL dynamics could be responsible and may impact individuals in different ways. In humans many lifestyle factors such as socio-economic status, psychological stress, smoking and body mass index are known to influence telomere length (Cherkas et al., 2008; Epel, Daubenmier, Moskowitz, Folkman, & Blackburn, 2009; Huzen et al., 2014; Kim et al., 2009; Shalev et al., 2012). In animal studies chronic infection and exposure to stressful situations have been linked to telomere shortening (Asghar, Hasselquist, Zehtindjiev, Westerdahl, & Bensch, 2015; Nettle et al., 2013; Reichert et al., 2014). Whilst we were unable to identify variables other than age which predicted within-individual TL change in our study system, individuals with a greater propensity to lose TL over time had shorter productive lives, implying changes in TL reflect important environmental or physiological variation linked to health. From a telomere biology perspective, it would be interesting to investigate changes not only in mean TL but also in its variance over time, preferentially using a technique like STELA or TESLA that allows the detection of critically short telomeres as well (Baird, Rowson, Wynford-Thomas, & Kipling, 2003; Lai et al., 2017). However, if the focus is to find good markers for animal breeding and health those techniques are not economically practical.

Our results have important implications for the application of TL as a biomarker within animal breeding and health. Productive lifespan describes the time from birth to culling and differs from measures of “longevity” or “lifespan” in human or wild animal studies because the animals are culled for reasons relating mainly to health and productivity. From an animal production perspective, a biomarker which could predict the productive lifespan of livestock from a relatively early age would represent an exciting target phenotype for use in selective breeding if heritable or for use in animal welfare monitoring if influenced by the environment. In the present study we showed that RLTL attrition predicted productive lifespan. The fact that both early life and lifelong RLTL shortening predicts earlier culling does support recent suggestions that TL may be useful as a biomarker of stress and welfare in livestock (Melissa Bateson, 2016; Brown, Dechow, Liu, Harvatine, & Ott, 2012). Many previous studies in birds and humans have found that experience of past stress is associated with shorter telomeres (Angelier, Vleck, Holberton, & Marra, 2013; Asghar et al., 2015; Cherkas et al., 2006; Entringer et al., 2011; Entringer et al., 2013; Epel et al., 2004; Haussmann, Longenecker, Marchetto, Juliano, & Bowden, 2012; Nettle et al., 2017). In our study, cattle showing RLTL attrition may be those that have experienced some form of recent immediate stress or minor infection. Repeated bouts of such stress, reflected by higher average RLTL attrition through life, may result in the development of health issues that lead to culling. Whilst the results of the present study suggest that longitudinal RLTL change itself may be a useful biomarker of health and welfare in livestock, a key area for further study is the causes of variation in RLTL change, as these may help us to understand sources of physiological stress for domestic animals and help improve animal welfare and productivity.

Acknowledgements

We thank everyone involved in collecting blood samples on the Crichton farm, particularly Hannah Dykes and Isla McCubbin. We also thank the SynthSys facility at the University of Edinburgh, particularly Eliane Salvo-Chirnside, for their support with the liquid handling robot. Our work was funded by the Biotechnology and Biological Sciences Research Council (BB/L007312/1; http://www.bbsrc.ac.uk/) awarded to GB. LS received a studentship from Scotland’s Rural College (http://www.sruc.ac.uk/).

References

Angelier, F., Costantini, D., Blévin, P., & Chastel, O. (2017). Do glucocorticoids mediate the link between environmental conditions and telomere dynamics in wild vertebrates? A review. General and Comparative Endocrinology, 256, 99–111. https://doi.org/10.1016/j.ygcen.2017.07.007

Angelier, F., Vleck, C. M., Holberton, R. L., & Marra, P. P. (2013). Telomere length, non-breeding habitat and return rate in male American redstarts. Functional Ecology, 27(2), 342–350. https://doi.org/10.1111/1365-2435.12041

Asghar, M, Hasselquist, D., Zehtindjiev, P., Westerdahl, H., & Bensch, S. (2015). Hidden costs of infection: Chronic malaria accelerates telomere degradation and senescence in wild birds, 347(6220), 9–12.

Asghar, Muhammad, Bensch, S., Tarka, M., Hansson, B., Hasselquist, D., & Asghar, M. (2015). Maternal and genetic factors determine early life telomere length. Proceedings of the Royal Society B.

Atema, E., Mulder, E., Dugdale, H. L., Briga, M., van Noordwijk, A. J., & Verhulst, S. (2015). Heritability of telomere length in the Zebra Finch. Journal of Ornithology, 156(4), 1113–1123. https://doi.org/10.1007/s10336-015-1212-7

Aubert, G., & Lansdorp, P. M. (2008). Telomeres and aging. Physiological Reviews, 88(2), 557–579. https://doi.org/10.1152/physrev.00026.2007

Aviv, A., & Shay, J. W. (2018). Reflections on telomere dynamics and ageing-related diseases in humans. Philosophical Transactions of the Royal Society of London. Series B, Biological Sciences, 373(1741), 20160436. https://doi.org/10.1098/rstb.2016.0436

Baerlocher, G. M., Rice, K., Vulto, I., & Lansdorp, P. M. (2007). Longitudinal data on telomere length in leukocytes from newborn baboons support a marked drop in stem cell turnover around 1 year of age. Aging Cell, 6(1), 121–123. https://doi.org/10.1111/j.1474-9726.2006.00254.x

Baird, D. M., Rowson, J., Wynford-Thomas, D., & Kipling, D. (2003). Extensive allelic variation and ultrashort telomeres in senescent human cells. Nature Genetics, 33(2), 203–207. https://doi.org/10.1038/ng1084

Bakaysa, S. L., Mucci, L. a., Slagboom, P. E., Boomsma, D. I., Mcclearn, G. E., Johansson, B., & Pedersen, N. L. (2007). Telomere length predicts survival independent of genetic influences. Aging Cell, 6(6), 769–774. https://doi.org/10.1111/j.1474-9726.2007.00340.x

Barrett, E. L. B., Burke, T. a., Hammers, M., Komdeur, J., & Richardson, D. S. (2013). Telomere length and dynamics predict mortality in a wild longitudinal study. Molecular Ecology, 22(1), 249–259. https://doi.org/10.1111/mec.12110

Bates, D., Mächler, M., Bolker, B., & Walker, S. (2015). Fitting Linear Mixed-Effects Models Using lme4. Journal of Statistical Software, 67(1). https://doi.org/10.18637/jss.v067.i01

Bateson, M, & Nettle, D. (2016). The telomere lengthening conundrum - it could be biology. Aging Cell, 2016, 1–8. https://doi.org/10.1111/acel.12555

Bateson, Melissa. (2016). Cumulative stress in research animals: Telomere attrition as a biomarker in a welfare context? BioEssays, 38(2), 201–212. https://doi.org/10.1002/bies.201500127

Benetos, A., Kark, J. D., Susser, E., Kimura, M., Sinnreich, R., Chen, W., … Aviv, A. (2013). Tracking and fixed ranking of leukocyte telomere length across the adult life course. Aging Cell, 12(4), 615–621. https://doi.org/10.1111/acel.12086

Bichet, C., Bouwhuis, S., Bauch, C., Verhulst, S., Becker, P. H., & Vedder, O. (2019). Telomere length is repeatable, shortens with age and reproductive success, and predicts remaining lifespan in a long‐lived seabird. Molecular Ecology, 2020(29), 429–441. https://doi.org/https://doi.org/10.1111/mec.1533

Bize, P., Criscuolo, F., Metcalfe, N. B., Nasir, L., & Monaghan, P. (2009). Telomere dynamics rather than age predict life expectancy in the wild. Proceedings of the Royal Society of London B, 276(1662), 1679–1683. https://doi.org/10.1098/rspb.2008.1817

Blackburn, E., & Gall, J. (1978). A Tandemly Repeated Sequence at the Termini of Extrachromosomal Ribbosomal RNA Genes in Tetrahymena. Journal of Molecular Biology, (120), 33–53.

Boonekamp, J., Mulder, G., Salomons, H. M., Dijkstra, C., & Verhulst, S. (2014). Nestling telomere shortening, but not telomere length, reflects developmental stress and predicts survival in wild birds. Proceedings of the Royal Society B: Biological Sciences, 281(April), 20133287. https://doi.org/http://dx.doi.org/10.1098/rspb.2013.3287

Boonekamp, J., Simons, M. J. P., Hemerik, L., & Verhulst, S. (2013). Telomere length behaves as biomarker of somatic redundancy rather than biological age. Aging Cell, 12(2), 330–332. https://doi.org/10.1111/acel.12050

Boonekamp, J., Mulder, G., Salomons, H., Dijkstra, C., & Verhulst, S. (2014). Nestling telomere shortening, but not telomere length, reflects developmental stress and predicts survival in wild birds. Proceedings of the Royal Society B: Biological Sciences, 281(April), 20133287. https://doi.org/http://dx.doi.org/10.1098/rspb.2013.3287

Broer, L., Codd, V., Nyholt, D. R., Deelen, J., Mangino, M., Willemsen, G., … Boomsma, D. I. (2013). Meta-analysis of telomere length in 19 713 subjects reveals high heritability, stronger maternal inheritance and a paternal age effect. European Journal of Human Genetics, 21(10), 1163–1168. https://doi.org/10.1038/ejhg.2012.303

Brown, D. E., Dechow, C. D., Liu, W. S., Harvatine, K. J., & Ott, T. L. (2012). Hot topic: Association of telomere length with age, herd, and culling in lactating Holsteins. Journal of Dairy Science, 95(11), 6384–6387. https://doi.org/10.3168/jds.2012-5593

Cawthon, R., Smith, K., O’Brien, E., Sivatchenko, A., & Kerber, R. (2003). Association between telomere length in blood and mortality in people aged 60 years or older. The Lancet, 361(9355), 393–395. https://doi.org/10.1016/S0140-6736(03)12384-7

Chen, W., Kimura, M., Kim, S., Cao, X., Srinivasan, S. R., Berenson, G. S., … Aviv, A. (2011). Longitudinal versus cross-sectional evaluations of leukocyte telomere length dynamics: Age-dependent telomere shortening is the rule. Journals of Gerontology - Series A Biological Sciences and Medical Sciences, 66 A(3), 312–319. https://doi.org/10.1093/gerona/glq223

Cherkas, L F, Aviv, a, Valdes, a M., Hunkin, J. L., Gardner, J. P., Surdulescu, G. L., … Spector, T. D. (2006). The effects of social status on biological aging as measured by white-blood-cell telomere length. Aging Cell, 5(5), 361–365. https://doi.org/10.1111/j.1474-9726.2006.00222.x

Cherkas, Lynn F, Hunkin, J. L., Kato, B. S., Richards, J. B., Gardner, J. P., Surdulescu, G. L., … Aviv, A. (2008). The association between physical activity in leisure time and leukocyte telomere length. Archives of Internal Medicine, 168(2), 154–158. https://doi.org/10.1001/archinternmed.2007.39

De Lange, T. (2005). Shelterin: The protein complex that shapes and safeguards human telomeres. Genes and Development, 19(18), 2100–2110. https://doi.org/10.1101/gad.1346005

Dugdale, H. L., & Richardson, D. S. (2018). Heritability of telomere variation: it’s all about the environment! Philosophical Transactions of the Royal Society B, 373(1741), 20160450. https://doi.org/10.1098/rstb.2016.0450

Entringer, S., Epel, E. S., Kumsta, R., Lin, J., Hellhammer, D. H., Blackburn, E. H., … Wadhwa, P. D. (2011). Stress exposure in intrauterine life is associated with shorter telomere length in young adulthood. Proceedings of the National Academy of Sciences, 108(33), E513–E518. https://doi.org/10.1073/pnas.1107759108

Entringer, Sonja, Epel, E. S., Lin, J., Buss, C., Shahbaba, B., Blackburn, E. H., … Wadhwa, P. D. (2013). Maternal psychosocial stress during pregnancy is associated with newborn leukocyte telomere length. American Journal of Obstetrics and Gynecology, 208(2), 134.e1-134.e7. https://doi.org/10.1016/j.ajog.2012.11.033

Epel, E., Daubenmier, J., Moskowitz, J. T., Folkman, S., & Blackburn, E. (2009). Can meditation slow rate of cellular aging? Cognitive stress, mindfulness, and telomeres. Annals of the New York Academy of Sciences, 1172, 34–53. https://doi.org/10.1111/j.1749-6632.2009.04414.x

Epel, E. S., Blackburn, E. H., Lin, J., Dhabhar, F. S., Adler, N. E., Morrow, J. D., & Cawthon, R. M. (2004). Accelerated telomere shortening in response to life stress. Proceedings of the National Academy of Sciences, 101(49), 17312–17315. https://doi.org/10.1073/pnas.0407162101

Fairlie, J., Holland, R., Pilkington, J. G., Pemberton, J. M., Harrington, L., & Nussey, D. H. (2015). Lifelong leukocyte telomere dynamics and survival in a free-living mammal. Aging Cell, 140–148. https://doi.org/10.1111/acel.12417

Farzaneh-Far, R., Lin, J., Epel, E., Lapham, K., Blackburn, E., & Whooley, M. A. (2010). Telomere length trajectory and its determinants in persons with coronary artery disease: Longitudinal findings from the heart and soul study. PLoS ONE, 5(1). https://doi.org/10.1371/journal.pone.0008612

Gardner, J. P., Li, S., Srinivasan, S. R., Chen, W., Kimura, M., Lu, X., … Aviv, A. (2005). Rise in insulin resistance is associated with escalated telomere attrition. Circulation, 111(17), 2171–2177. https://doi.org/10.1161/01.CIR.0000163550.70487.0B

Goglin, S. E., Farzaneh-Far, R., Epel, E. S., Lin, J., Blackburn, E. H., & Whooley, M. A. (2016). Change in leukocyte telomere length predicts mortality in patients with stable coronary heart disease from the heart and soul study. PLoS ONE, 11(10), 1–12. https://doi.org/10.1371/journal.pone.0160748

Gomes, N. M. V, Ryder, O. a., Houck, M. L., Charter, S. J., Walker, W., Forsyth, N. R., … Wright, W. E. (2011). Comparative biology of mammalian telomeres: Hypotheses on ancestral states and the roles of telomeres in longevity determination. Aging Cell, 10(5), 761–768. https://doi.org/10.1111/j.1474-9726.2011.00718.x

Harley, C., Futcher, A., & Greider, C. (1990). Telomeres shorten during ageing of human fibroblasts. Nature, 345(6274), 458–460. https://doi.org/10.1038/345458a0

Harrington, L., & Pucci, F. (2018). In medio stat virtus: unanticipated consequences of telomere dysequilibrium. Philosophical Transactions of the Royal Society B: Biological Sciences, 373(1741), 20160444. https://doi.org/10.1098/rstb.2016.0444

Haussmann, M. F., Longenecker, a. S., Marchetto, N. M., Juliano, S. a., & Bowden, R. M. (2012). Embryonic exposure to corticosterone modifies the juvenile stress response, oxidative stress and telomere length. Proceedings of the Royal Society B: Biological Sciences, 279(1732), 1447–1456. https://doi.org/10.1098/rspb.2011.1913

Haussmann, Mark F, Winkler, D. W., & Vleck, C. M. (2005). Longer telomeres associated with higher survival in birds. Biology Letters, 1(2), 212–214. https://doi.org/10.1098/rsbl.2005.0301

Heidinger, B. J., Blount, J. D., Boner, W., Griffiths, K., Metcalfe, N. B., & Monaghan, P. (2012). Telomere length in early life predicts lifespan. Proceedings of the National Academy of Sciences, 109(5), 1743–1748. https://doi.org/10.1073/pnas.1113306109

Huzen, J., Wong, L. S. M., van Veldhuisen, D. J., Samani, N. J., Zwinderman, A. H., Codd, V., … Van der Harst, P. (2014). Telomere length loss due to smoking and metabolic traits. Journal of Internal Medicine, 275(2), 155–163. https://doi.org/10.1111/joim.12149

Kark, J., Goldberger, N., Kimura, M., Sinnreich, R., & Aviv, A. (2012). Energy intake and leukocyte telomere length in young adults. The American Journal of Clinical Nutrition, (95), 479–487. https://doi.org/10.3945/ajcn.111.024521.1

Kim, S., Parks, C. G., DeRoo, L. A., Chen, H., Taylor, J. A., Cawthon, R. M., & Sandler, D. P. (2009). Obesity and weight gain in adulthood and telomere length. Cancer Epidemiology Biomarkers and Prevention, 18(3), 816–820. https://doi.org/10.1158/1055-9965.EPI-08-0935

Lai, T.-P., Zhang, N., Noh, J., Mender, I., Tedone, E., Huang, E., … Shay, J. W. (2017). A method for measuring the distribution of the shortest telomeres in cells and tissues. Nature Communications, 8(1), 1356. https://doi.org/10.1038/s41467-017-01291-z

Muñoz-lorente, M. A., Cano-martin, A. C., & Blasco, M. A. (2019). Mice with hyper-long telomeres show less metabolic aging and longer lifespans. Nature Communications, 10(4723), 1–14. https://doi.org/10.1038/s41467-019-12664-x

Nettle, D., Andrews, C., Reichert, S., Bedford, T., Kolenda, C., Parker, C., … Bateson, M. (2017). Early-life adversity accelerates cellular ageing and affects adult inflammation: Experimental evidence from the European starling. Scientific Reports, 7: 40794(December 2016), 1–10. https://doi.org/10.1038/srep40794

Nettle, D., Monaghan, P., Boner, W., Gillespie, R., & Bateson, M. (2013). Bottom of the heap: Having heavier competitors accelerates early-life telomere loss in the European starling, Sturnus vulgaris. PLoS ONE, 8(12). https://doi.org/10.1371/journal.pone.0083617

Nettle, D., Seeker, L., Nussey, D., Froy, H., & Bateson, M. (2019). Consequences of measurement error in qPCR telomere data: A simulation study. PLOS ONE, 14(5), e0216118. https://doi.org/10.1371/journal.pone.0216118

Nordfjäll, K., Svenson, U., Norrback, K. F., Adolfsson, R., Lenner, P., & Roos, G. (2009). The individual blood cell telomere attrition rate is telomere length dependent. PLoS Genetics, 5(2), 2–7. https://doi.org/10.1371/journal.pgen.1000375

Olovnikov, A. M. (1973). A theory of marginotomy. The incomplete copying of template margin in enzymic synthesis of polynucleotides and biological significance of the phenomenon. Journal of Theoretical Biology, 41(1), 181–190. https://doi.org/10.1016/0022-5193(73)90198-7

R Core Team. (2014). R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing. Retrieved from http://www.r-project.org/

Reichert, S., Stier, A., Zahn, S., Arrivé, M., Bize, P., Massemin, S., & Criscuolo, F. (2014). Increased brood size leads to persistent eroded telomeres. Frontiers in Ecology and Evolution, 2(April), 1–11. https://doi.org/10.3389/fevo.2014.00009

Seeker, L. A., Holland, R., Underwood, S., Fairlie, J., Psifidi, A., Ilska, J., … Nussey, D. H. (2016). Method specific calibration corrects for DNA extraction method effects on relative telomere length measurements by quantitative PCR. PLoS ONE, 11(10), 1–15. https://doi.org/10.1371/journal.pone.0164046

Seeker, L. A., Ilska, J., Psifidi, A., Wilbourn, R. V., Underwood, S. L., Fairlie, J., … Banos, G. (2018). Bovine telomere dynamics and the association between telomere length and productive lifespan. Scientific Reports.

Seeker, L. A., Ilska, J., Psifidi, A., Wilbourn, R. V, Underwood, L., Fairlie, J., … Banos, G. (2018). Longitudinal changes in telomere length and associated genetic parameters in dairy cattle analysed using random regression models. Plos One, 1–15. https://doi.org/10.1371/journal.pone.0192864

Seeker, L., Holland, B., Psifidi, A., Banos, B., & Nussey, D. (2015). A robust assay for the measurement of telomere length in dairy cattle. In BSAS annual conference (p. 066).

Shalev, I., Moffitt, T. E., Sugden, K., Williams, B., Houts, R. M., Danese, A., … Caspi, A. (2012). Exposure to violence during childhood is associated with telomere erosion from 5 to 10 years of age: a longitudinal study. Mol.Psychiatry, 18(5), 576–581. https://doi.org/10.1038/mp.2012.32

Spurgin, L. G., Bebbington, K., Fairfield, E. A., Hammers, M., Komdeur, J., Burke, T., … Richardson, D. S. (2018). Spatio-temporal variation in lifelong telomere dynamics in a long-term ecological study. Journal of Animal Ecology, 87(1), 187–198. https://doi.org/10.1111/1365-2656.12741

Steenstrup, T., Hjelmborg, J. V. B., Mortensen, L. H., Kimura, M., Christensen, K., & Aviv, A. (2013). Leukocyte telomere dynamics in the elderly. European Journal of Epidemiology, 28(2), 181–187. https://doi.org/10.1007/s10654-013-9780-4

Svenson, U., Nordfjäll, K., Baird, D., Roger, L., Osterman, P., Hellenius, M.-L., & Roos, G. (2011). Blood Cell Telomere Length Is a Dynamic Feature. PLOS ONE, 6(6), e21485. Retrieved from https://doi.org/10.1371/journal.pone.0021485

Therneau, T. M. (2015). A Package for Survival Analysis in S. Retrieved from https://cran.r-project.org/package=survival

Veerkamp, R. ., Simm, G., & Oldham, J. . (1994). Effects of interaction between genotype and feeding system on milk production, feed intake, efficiency and body tissue mobilization in dairy cows. Livestock Production Science, 39(3), 229–241. https://doi.org/10.1016/0301-6226(94)90202-X

Voillemot, M., Hine, K., Zahn, S., Criscuolo, F., Gustafsson, L., Doligez, B., & Bize, P. (2012). Effects of brood size manipulation and common origin on phenotype and telomere length in nestling collared flycatchers. BMC Ecology, 12(1), 17. https://doi.org/10.1186/1472-6785-12-17

Watson, J. D. (1972). Origin of concatemeric T7 DNA. Nature: New Biology, 239(94), 197–201. Retrieved from http://www.ncbi.nlm.nih.gov/pubmed/4507727

Whittemore, K., Vera, E., Martínez-nevado, E., Sanpera, C., & Blasco, M. A. (2019). Telomere shortening rate predicts species life span. Proceedings of the National Academy of Sciences, 116(30), 15122–15127. https://doi.org/10.1073/pnas.1902452116

Wickham, H. (2009). ggplot2: Elegant Graphics for Data Analysis. Springer-Verlag New York. Retrieved from http://ggplot2.org

Wilbourn, R. V., Moatt, J. P., Froy, H., Walling, C. A., Nussey, D. H., & Boonekamp, J. J. (2018). The relationship between telomere length and mortality risk in non-model vertebrate systems: a meta-analysis. Phil. Trans. R. Soc. B, 373(1741), 20160447. https://doi.org/10.1098/RSTB.2016.0447

Yuan, X., Kronström, M., Hellenius, M. L., Cederholm, T., & Xu, D. (2018). Longitudinal changes in leukocyte telomere length and mortality in elderly Swedish men. Aging, 10(10), 3005–3016.

Data Accessibility

The complete dataset will be accessible after acceptance on GitHub.

Author contributions

All authors designed the experiment together. LS, SU, RW, JF, RH and AB collected the blood samples and processed them in the lab. LS analysed the data with help from HF, JI, AP, GB and DN. LS wrote the initial draft of the manuscript, DN and GB edited. All authors were involved in revising the manuscript. BW, MC, GB and DN obtained funding for the project.

PART 2: THE ANALYSIS

Load libraries

library(ggplot2)
library(nlme)
library(lme4)
## Loading required package: Matrix
## 
## Attaching package: 'lme4'
## The following object is masked from 'package:nlme':
## 
##     lmList
library(lmerTest)
## 
## Attaching package: 'lmerTest'
## The following object is masked from 'package:lme4':
## 
##     lmer
## The following object is masked from 'package:stats':
## 
##     step
library(grid)
library(gridExtra)
library(mgcv) # for GAM
## This is mgcv 1.8-31. For overview type 'help("mgcv-package")'.
library(data.table)
library(stringr)

library(magrittr)
library(ggthemes)
library(plyr)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:plyr':
## 
##     arrange, count, desc, failwith, id, mutate, rename, summarise,
##     summarize
## The following objects are masked from 'package:data.table':
## 
##     between, first, last
## The following object is masked from 'package:gridExtra':
## 
##     combine
## The following object is masked from 'package:nlme':
## 
##     collapse
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(survival)
library(survminer)
## Loading required package: ggpubr
## Registered S3 methods overwritten by 'car':
##   method                          from
##   influence.merMod                lme4
##   cooks.distance.influence.merMod lme4
##   dfbeta.influence.merMod         lme4
##   dfbetas.influence.merMod        lme4
## 
## Attaching package: 'ggpubr'
## The following object is masked from 'package:plyr':
## 
##     mutate
library(ggsci)
#library(AICcmodavg)

Build colour palettes

# Make colour palette
mypal <- pal_npg("nrc", alpha = 0.7)(10)
mypal2 <-pal_tron("legacy", alpha = 0.7)(7)
mypal3<-pal_lancet("lanonc", alpha = 0.7)(9)
mypal4<-pal_simpsons(palette = c("springfield"), alpha = 0.7)(16)

mycoloursP<- c(mypal, mypal2, mypal3, mypal4)

Read in data

data<- read.csv("/Users/luise/Documents/Work/CattleTelomereLength/LongTelomereChangePaper/git/TelomereChangeInDairyCattle/LongTelomereDataCattle.csv")


# create data which only lists each animal ID once (needed for some counts/ histograms below)
uData<- data[!duplicated(data$recoded_id),]

# create data that excludes all first measurements without t-1 time point. This is used for investigations of change
redData<-subset(data, data$residual_rltl_tminus1 != "NULL")

init_change_dat<-read.csv("/Users/luise/Documents/Work/CattleTelomereLength/LongTelomereChangePaper/git/TelomereChangeInDairyCattle/first_year_long_tel_cattle.csv")


coxData<-read.csv("/Users/luise/Documents/Work/CattleTelomereLength/LongTelomereChangePaper/git/TelomereChangeInDairyCattle/cox_data.csv")

coxDataR<-subset(coxData, coxData$amp_305_lactation != "NULL")

Fig. 1

cor.test(redData$residual_rltl_t,as.numeric(paste(redData$residual_rltl_tminus1)))
## 
##  Pearson's product-moment correlation
## 
## data:  redData$residual_rltl_t and as.numeric(paste(redData$residual_rltl_tminus1))
## t = 13.058, df = 1018, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.3249385 0.4301473
## sample estimates:
##       cor 
## 0.3787659
SubsPlot<-ggplot(redData, aes(x=as.numeric(paste(residual_rltl_tminus1)), y=redData$residual_rltl_t))+geom_point()+xlab("RLTL at time t-1")+ylab("RLTL at time t")+theme_classic(20)+geom_abline(intercept=0, slope=1, colour="red")+ stat_smooth(method=lm, formula=y~x)+
  annotate("text", x = -0.4, y = 0.5, label = "(B)", fontface =2, size = 10)


agePlot<- ggplot(redData, aes(x=as.factor(redData$age_y), y=as.numeric(paste(redData$diff_rltl_res)), colour = as.factor(redData$age_y)))+geom_jitter()+theme_classic(20)+ xlab("Age in years")+ ylab("RLTL change")+ geom_boxplot()+
  annotate("text", x = 1, y = 0.5, label = "(A)", fontface =2, size = 10) + scale_color_manual(values=mycoloursP[16:40])+ theme(legend.position="none")

grid.arrange(agePlot, SubsPlot, ncol=2)
## Warning: Use of `redData$age_y` is discouraged. Use `age_y` instead.
## Warning: Use of `redData$diff_rltl_res` is discouraged. Use `diff_rltl_res`
## instead.
## Warning: Use of `redData$age_y` is discouraged. Use `age_y` instead.

## Warning: Use of `redData$age_y` is discouraged. Use `age_y` instead.
## Warning: Use of `redData$diff_rltl_res` is discouraged. Use `diff_rltl_res`
## instead.
## Warning: Use of `redData$age_y` is discouraged. Use `age_y` instead.
## Warning: Use of `redData$residual_rltl_t` is discouraged. Use `residual_rltl_t`
## instead.

## Warning: Use of `redData$residual_rltl_t` is discouraged. Use `residual_rltl_t`
## instead.

Fig. 2

surv_plotlist<-list()
surv_plotlist[[1]] <- init_ch_plot 
surv_plotlist[[2]] <- mean_abs_rltl_change
surv_plotlist[[3]] <- meanRLTL
surv_plotlist[[4]] <- mean_rltl_change



arrange_ggsurvplots(surv_plotlist, print = TRUE, ncol = 2, nrow = 2)

calculate sample number per animal

count<- nrow(data)
unique_count<- length(unique(data$recoded_id))
sample_year_summ<-summary(as.factor(data$sample_year))

print(paste("There are ", count, " relative leukocyte telomere length (RLTL) measurements of ", unique_count, " animals included in the study.", sep =""))
## [1] "There are 1325 relative leukocyte telomere length (RLTL) measurements of 305 animals included in the study."

Distibution of sample years

sample_year_summ
## 2008 2009 2010 2011 2012 2013 2014 
##   22   91  179  247  266  272  248

Plot sample number per animal

hist1<-ggplot(uData, aes(x = uData$SamplesPerAnimal)) + 
    geom_histogram(colour= "black", fill = mycoloursP[7], binwidth = 1) + xlab("Number of samples per animal") +  
  geom_vline(xintercept= mean(as.numeric(paste(uData$SamplesPerAnimal))), colour="red")+ 
  annotate("text", x = 1, y = 100, label = "(A)", fontface =2, size = 8) + theme_classic(20)+
  scale_x_continuous(breaks=c(1:8))

#hist1

Plot distribution of telomere residuals

hist2<-ggplot(data, aes(x = data$residual_rltl_t)) + 
    geom_histogram(colour= "black", fill = mycoloursP[12]) + xlab("Adjusted RLTL") + 
  geom_vline(xintercept= mean(as.numeric(paste(data$residual_rltl_t))), colour="red") + 
  annotate("text", x = -0.4, y = 120, label = "(B)", fontface =2, size = 8)+
  theme_classic(20)+
    stat_function(
        fun = function(x, mean, sd, n){
             n/25* dnorm(x = x, mean = mean, sd = sd)
        }, 
        args = with(data, c(mean = mean(data$residual_rltl_t), sd = sd(data$residual_rltl_t), n
= length(data$residual_rltl_t)))
    ) 

#hist2
# statistical tests for normal distribution of RTL measurements
shapiro.test(data$residual_rltl_t)
## 
##  Shapiro-Wilk normality test
## 
## data:  data$residual_rltl_t
## W = 0.989, p-value = 1.967e-08
ks.test(data$residual_rltl_t, pnorm)
## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  data$residual_rltl_t
## D = 0.36744, p-value < 2.2e-16
## alternative hypothesis: two-sided
qqnorm(data$residual_rltl_t)

Plot difference between RLTL res at time t and time t-1

diff<- as.numeric(paste(redData$diff_rltl_res))

hist3<-ggplot(redData, aes(x = as.numeric(paste(redData$diff_rltl_res)))) + 
    geom_histogram(colour= "black", fill = mycoloursP[13]) + xlab("Adjusted RLTL change") +  
  geom_vline(xintercept= mean(as.numeric(paste(redData$diff_rltl_res))), colour="red") +
  annotate("text", x = -0.6, y = 95, label = "(C)", fontface =2, size= 8)+
  theme_classic(20)+
    stat_function(
        fun = function(x, mean, sd, n){
             n/23* dnorm(x = x, mean = mean, sd = sd)
        }, 
        args = with(redData, c(mean = mean(as.numeric(paste(redData$diff_rltl_res))), sd = sd(as.numeric(paste(redData$diff_rltl_res))), n
= length(as.numeric(paste(redData$diff_rltl_res))))))

#hist3
grid.arrange(hist1, hist2, hist3, ncol = 3)
## Warning: Use of `uData$SamplesPerAnimal` is discouraged. Use `SamplesPerAnimal`
## instead.
## Warning: Use of `data$residual_rltl_t` is discouraged. Use `residual_rltl_t`
## instead.

## Warning: Use of `data$residual_rltl_t` is discouraged. Use `residual_rltl_t`
## instead.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning: Use of `redData$diff_rltl_res` is discouraged. Use `diff_rltl_res`
## instead.
## Warning: Use of `redData$diff_rltl_res` is discouraged. Use `diff_rltl_res`
## instead.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

# statistical tests for normal distribution of difference in RTL residuals
shapiro.test(diff)
## 
##  Shapiro-Wilk normality test
## 
## data:  diff
## W = 0.998, p-value = 0.2655
ks.test(diff, pnorm)
## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  diff
## D = 0.35414, p-value < 2.2e-16
## alternative hypothesis: two-sided
qqnorm(diff)

# test if mean RLTL change is significantly different from zero:

t.test(as.numeric(paste(redData$diff_rltl_res)), mu=0, alternative="less", conf.level=0.99)
## 
##  One Sample t-test
## 
## data:  as.numeric(paste(redData$diff_rltl_res))
## t = -5.8438, df = 1019, p-value = 3.43e-09
## alternative hypothesis: true mean is less than 0
## 99 percent confidence interval:
##         -Inf -0.01949685
## sample estimates:
##   mean of x 
## -0.03242526

Shortening and elongating:

shortening<-subset(redData, as.numeric(paste(redData$diff_rltl_res)) < 0)
elongating<-subset(redData, as.numeric(paste(redData$diff_rltl_res))>0)

percent_shortening<-nrow(shortening)/nrow(redData)*100
percent_shortening
## [1] 56.47059
percent_elongating<-nrow(elongating)/nrow(redData)*100
percent_elongating
## [1] 43.52941

Initial full model of RLTL change

init_mod<-lmer(as.numeric(paste(redData$diff_rltl_res))~ redData$age_y + redData$feed_group + redData$genetic_group+ as.factor(redData$birth_year)+ as.factor(redData$sample_year) +  as.numeric(paste(redData$sample_interval))+ redData$health_event_2week +(1|redData$recoded_id), na.action=na.exclude)
## boundary (singular) fit: see ?isSingular
summary(init_mod)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: 
## as.numeric(paste(redData$diff_rltl_res)) ~ redData$age_y + redData$feed_group +  
##     redData$genetic_group + as.factor(redData$birth_year) + as.factor(redData$sample_year) +  
##     as.numeric(paste(redData$sample_interval)) + redData$health_event_2week +  
##     (1 | redData$recoded_id)
## 
## REML criterion at convergence: -599.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.4980 -0.6827  0.0248  0.6712  3.5498 
## 
## Random effects:
##  Groups             Name        Variance Std.Dev.
##  redData$recoded_id (Intercept) 0.00000  0.0000  
##  Residual                       0.02953  0.1718  
## Number of obs: 1020, groups:  redData$recoded_id, 305
## 
## Fixed effects:
##                                              Estimate Std. Error         df
## (Intercept)                                 2.376e-01  1.737e-01  1.002e+03
## redData$age_y                               2.219e-02  1.385e-02  1.002e+03
## redData$feed_group2                        -8.381e-03  1.223e-02  1.002e+03
## redData$feed_groupNULL                     -6.651e-03  1.627e-02  1.002e+03
## redData$genetic_groupS                     -1.025e-02  1.144e-02  1.002e+03
## as.factor(redData$birth_year)2009          -1.431e-02  2.062e-02  1.002e+03
## as.factor(redData$birth_year)2010          -1.231e-02  2.855e-02  1.002e+03
## as.factor(redData$birth_year)2011          -6.343e-03  3.863e-02  1.002e+03
## as.factor(redData$birth_year)2012          -2.192e-02  5.533e-02  1.002e+03
## as.factor(redData$birth_year)2013          -1.466e-01  1.114e-01  1.002e+03
## as.factor(redData$birth_year)2014          -2.415e-01  1.925e-01  1.002e+03
## as.factor(redData$sample_year)2010         -3.639e-01  1.748e-01  1.002e+03
## as.factor(redData$sample_year)2011         -3.223e-01  1.763e-01  1.002e+03
## as.factor(redData$sample_year)2012         -3.008e-01  1.784e-01  1.002e+03
## as.factor(redData$sample_year)2013         -2.767e-01  1.824e-01  1.002e+03
## as.factor(redData$sample_year)2014         -3.003e-01  1.886e-01  1.002e+03
## as.numeric(paste(redData$sample_interval)) -2.685e-06  5.434e-05  1.002e+03
## redData$health_event_2week                  1.764e-02  2.447e-02  1.002e+03
##                                            t value Pr(>|t|)  
## (Intercept)                                  1.368   0.1715  
## redData$age_y                                1.602   0.1095  
## redData$feed_group2                         -0.685   0.4933  
## redData$feed_groupNULL                      -0.409   0.6828  
## redData$genetic_groupS                      -0.896   0.3706  
## as.factor(redData$birth_year)2009           -0.694   0.4880  
## as.factor(redData$birth_year)2010           -0.431   0.6663  
## as.factor(redData$birth_year)2011           -0.164   0.8696  
## as.factor(redData$birth_year)2012           -0.396   0.6920  
## as.factor(redData$birth_year)2013           -1.316   0.1886  
## as.factor(redData$birth_year)2014           -1.255   0.2099  
## as.factor(redData$sample_year)2010          -2.082   0.0376 *
## as.factor(redData$sample_year)2011          -1.829   0.0677 .
## as.factor(redData$sample_year)2012          -1.686   0.0921 .
## as.factor(redData$sample_year)2013          -1.517   0.1295  
## as.factor(redData$sample_year)2014          -1.592   0.1117  
## as.numeric(paste(redData$sample_interval))  -0.049   0.9606  
## redData$health_event_2week                   0.721   0.4711  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation matrix not shown by default, as p = 18 > 12.
## Use print(x, correlation=TRUE)  or
##     vcov(x)        if you need it
## convergence code: 0
## boundary (singular) fit: see ?isSingular
anova(init_mod)
## Type III Analysis of Variance Table with Satterthwaite's method
##                                             Sum Sq  Mean Sq NumDF DenDF F value
## redData$age_y                              0.07576 0.075762     1  1002  2.5657
## redData$feed_group                         0.01465 0.007327     2  1002  0.2481
## redData$genetic_group                      0.02369 0.023694     1  1002  0.8024
## as.factor(redData$birth_year)              0.12877 0.021462     6  1002  0.7268
## as.factor(redData$sample_year)             0.42749 0.085498     5  1002  2.8954
## as.numeric(paste(redData$sample_interval)) 0.00007 0.000072     1  1002  0.0024
## redData$health_event_2week                 0.01535 0.015351     1  1002  0.5199
##                                             Pr(>F)  
## redData$age_y                              0.10952  
## redData$feed_group                         0.78031  
## redData$genetic_group                      0.37059  
## as.factor(redData$birth_year)              0.62806  
## as.factor(redData$sample_year)             0.01329 *
## as.numeric(paste(redData$sample_interval)) 0.96060  
## redData$health_event_2week                 0.47107  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
AIC(init_mod)
## [1] -559.2719
# There is a warning that this model is singular which is due to the random effect variance structure. It means that the variance due to the animal ID is not significantly different from 0. Although this may be a problem for some models, I think in this case it is not. It simply means that there is no variation between animals in the RLTL change and that is to be expected here, because the variance within animals is much larger than the variance between animals. 

# Sample year is statistically significant, but it might only be signifgicant due to the design of our study. Becasue we wanted to have repeat measurements for all animals with the first measurement taken with the first 15 days of life, sampling in the first year only included calves which have been shown before to have longer telomeres than adults and to shorten them in their first year of life. So is the year effect simply driven by that? If so, it should vanish after excluding samples taken within the two years

Full model of RLTL change excluding first two years

summary(as.factor(data$sample_year))
## 2008 2009 2010 2011 2012 2013 2014 
##   22   91  179  247  266  272  248
test_dat<-subset(data, data$sample_year != 2008)
test_dat<-subset(test_dat, test_dat$sample_year != 2009)

long_dat<-subset(test_dat, test_dat$diff_rltl_res != "NULL")

late_mod<-lmer(as.numeric(paste(long_dat$diff_rltl_res))~ long_dat$age_y + long_dat$feed_group + long_dat$genetic_group + as.factor(long_dat$birth_year)+ as.factor(long_dat$sample_year) +  as.numeric(paste(long_dat$sample_interval))+ long_dat$health_event_2week +(1|long_dat$recoded_id), na.action=na.exclude)
## boundary (singular) fit: see ?isSingular
summary(late_mod)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: as.numeric(paste(long_dat$diff_rltl_res)) ~ long_dat$age_y +  
##     long_dat$feed_group + long_dat$genetic_group + as.factor(long_dat$birth_year) +  
##     as.factor(long_dat$sample_year) + as.numeric(paste(long_dat$sample_interval)) +  
##     long_dat$health_event_2week + (1 | long_dat$recoded_id)
## 
## REML criterion at convergence: -599.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.4980 -0.6836  0.0252  0.6717  3.5498 
## 
## Random effects:
##  Groups              Name        Variance Std.Dev.
##  long_dat$recoded_id (Intercept) 0.00000  0.0000  
##  Residual                        0.02953  0.1718  
## Number of obs: 1019, groups:  long_dat$recoded_id, 305
## 
## Fixed effects:
##                                               Estimate Std. Error         df
## (Intercept)                                 -1.263e-01  3.879e-02  1.002e+03
## long_dat$age_y                               2.219e-02  1.385e-02  1.002e+03
## long_dat$feed_group2                        -8.381e-03  1.223e-02  1.002e+03
## long_dat$feed_groupNULL                     -6.651e-03  1.627e-02  1.002e+03
## long_dat$genetic_groupS                     -1.025e-02  1.144e-02  1.002e+03
## as.factor(long_dat$birth_year)2009          -1.431e-02  2.062e-02  1.002e+03
## as.factor(long_dat$birth_year)2010          -1.231e-02  2.855e-02  1.002e+03
## as.factor(long_dat$birth_year)2011          -6.343e-03  3.863e-02  1.002e+03
## as.factor(long_dat$birth_year)2012          -2.192e-02  5.533e-02  1.002e+03
## as.factor(long_dat$birth_year)2013          -1.466e-01  1.114e-01  1.002e+03
## as.factor(long_dat$birth_year)2014          -2.415e-01  1.925e-01  1.002e+03
## as.factor(long_dat$sample_year)2011          4.158e-02  2.644e-02  1.002e+03
## as.factor(long_dat$sample_year)2012          6.309e-02  3.329e-02  1.002e+03
## as.factor(long_dat$sample_year)2013          8.721e-02  4.463e-02  1.002e+03
## as.factor(long_dat$sample_year)2014          6.360e-02  6.047e-02  1.002e+03
## as.numeric(paste(long_dat$sample_interval)) -2.685e-06  5.434e-05  1.002e+03
## long_dat$health_event_2week                  1.764e-02  2.447e-02  1.002e+03
##                                             t value Pr(>|t|)   
## (Intercept)                                  -3.256  0.00117 **
## long_dat$age_y                                1.602  0.10952   
## long_dat$feed_group2                         -0.685  0.49328   
## long_dat$feed_groupNULL                      -0.409  0.68275   
## long_dat$genetic_groupS                      -0.896  0.37059   
## as.factor(long_dat$birth_year)2009           -0.694  0.48796   
## as.factor(long_dat$birth_year)2010           -0.431  0.66631   
## as.factor(long_dat$birth_year)2011           -0.164  0.86961   
## as.factor(long_dat$birth_year)2012           -0.396  0.69205   
## as.factor(long_dat$birth_year)2013           -1.316  0.18862   
## as.factor(long_dat$birth_year)2014           -1.255  0.20990   
## as.factor(long_dat$sample_year)2011           1.573  0.11605   
## as.factor(long_dat$sample_year)2012           1.895  0.05835 . 
## as.factor(long_dat$sample_year)2013           1.954  0.05097 . 
## as.factor(long_dat$sample_year)2014           1.052  0.29309   
## as.numeric(paste(long_dat$sample_interval))  -0.049  0.96060   
## long_dat$health_event_2week                   0.721  0.47107   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation matrix not shown by default, as p = 17 > 12.
## Use print(x, correlation=TRUE)  or
##     vcov(x)        if you need it
## convergence code: 0
## boundary (singular) fit: see ?isSingular
anova(late_mod)
## Type III Analysis of Variance Table with Satterthwaite's method
##                                               Sum Sq  Mean Sq NumDF DenDF
## long_dat$age_y                              0.075762 0.075762     1  1002
## long_dat$feed_group                         0.014654 0.007327     2  1002
## long_dat$genetic_group                      0.023694 0.023694     1  1002
## as.factor(long_dat$birth_year)              0.128774 0.021462     6  1002
## as.factor(long_dat$sample_year)             0.312340 0.078085     4  1002
## as.numeric(paste(long_dat$sample_interval)) 0.000072 0.000072     1  1002
## long_dat$health_event_2week                 0.015351 0.015351     1  1002
##                                             F value  Pr(>F)  
## long_dat$age_y                               2.5657 0.10952  
## long_dat$feed_group                          0.2481 0.78031  
## long_dat$genetic_group                       0.8024 0.37059  
## as.factor(long_dat$birth_year)               0.7268 0.62806  
## as.factor(long_dat$sample_year)              2.6443 0.03235 *
## as.numeric(paste(long_dat$sample_interval))  0.0024 0.96060  
## long_dat$health_event_2week                  0.5199 0.47107  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
AIC(late_mod)
## [1] -561.2719
# sample year remains statistically significant. 

Reduced model of RLTL change

red_mod<-lmer(as.numeric(paste(redData$diff_rltl_res))~ redData$age_y + redData$feed_group + redData$genetic_group+ as.factor(redData$sample_year) +  as.numeric(paste(redData$sample_interval)) + (1|redData$recoded_id), na.action=na.exclude)
## boundary (singular) fit: see ?isSingular
summary(red_mod)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: 
## as.numeric(paste(redData$diff_rltl_res)) ~ redData$age_y + redData$feed_group +  
##     redData$genetic_group + as.factor(redData$sample_year) +  
##     as.numeric(paste(redData$sample_interval)) + (1 | redData$recoded_id)
## 
## REML criterion at convergence: -627.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.5005 -0.6953  0.0197  0.6760  3.5199 
## 
## Random effects:
##  Groups             Name        Variance Std.Dev.
##  redData$recoded_id (Intercept) 0.00000  0.0000  
##  Residual                       0.02947  0.1717  
## Number of obs: 1020, groups:  redData$recoded_id, 305
## 
## Fixed effects:
##                                              Estimate Std. Error         df
## (Intercept)                                 2.278e-01  1.724e-01  1.009e+03
## redData$age_y                               2.571e-02  6.026e-03  1.009e+03
## redData$feed_group2                        -1.045e-02  1.206e-02  1.009e+03
## redData$feed_groupNULL                     -1.499e-02  1.546e-02  1.009e+03
## redData$genetic_groupS                     -1.262e-02  1.134e-02  1.009e+03
## as.factor(redData$sample_year)2010         -3.693e-01  1.738e-01  1.009e+03
## as.factor(redData$sample_year)2011         -3.298e-01  1.736e-01  1.009e+03
## as.factor(redData$sample_year)2012         -3.087e-01  1.735e-01  1.009e+03
## as.factor(redData$sample_year)2013         -2.882e-01  1.740e-01  1.009e+03
## as.factor(redData$sample_year)2014         -3.190e-01  1.750e-01  1.009e+03
## as.numeric(paste(redData$sample_interval))  1.029e-05  5.207e-05  1.009e+03
##                                            t value Pr(>|t|)    
## (Intercept)                                  1.321   0.1868    
## redData$age_y                                4.266 2.18e-05 ***
## redData$feed_group2                         -0.867   0.3864    
## redData$feed_groupNULL                      -0.969   0.3327    
## redData$genetic_groupS                      -1.113   0.2660    
## as.factor(redData$sample_year)2010          -2.125   0.0338 *  
## as.factor(redData$sample_year)2011          -1.900   0.0578 .  
## as.factor(redData$sample_year)2012          -1.779   0.0755 .  
## as.factor(redData$sample_year)2013          -1.657   0.0979 .  
## as.factor(redData$sample_year)2014          -1.823   0.0686 .  
## as.numeric(paste(redData$sample_interval))   0.198   0.8434    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) rdDt$_ rdD$_2 rD$_NU rdD$_S a.(D$_)2010 a.(D$_)2011
## redData$g_y  0.001                                                    
## rdDt$fd_gr2 -0.068 -0.099                                             
## rdDt$f_NULL -0.042  0.114  0.370                                      
## rdDt$gntc_S -0.064  0.090 -0.029  0.243                               
## a.(D$_)2010 -0.987 -0.034  0.032  0.010  0.021                        
## a.(D$_)2011 -0.988 -0.052  0.034  0.001  0.018  0.991                 
## a.(D$_)2012 -0.989 -0.067  0.037  0.006  0.019  0.992       0.995     
## a.(D$_)2013 -0.986 -0.094  0.038 -0.001  0.015  0.991       0.995     
## a.(D$_)2014 -0.980 -0.130  0.039 -0.005  0.009  0.988       0.992     
## as.n((D$_)) -0.008 -0.033  0.047 -0.034  0.071 -0.099      -0.105     
##             a.(D$_)2012 a.(D$_)2013 a.(D$_)2014
## redData$g_y                                    
## rdDt$fd_gr2                                    
## rdDt$f_NULL                                    
## rdDt$gntc_S                                    
## a.(D$_)2010                                    
## a.(D$_)2011                                    
## a.(D$_)2012                                    
## a.(D$_)2013  0.996                             
## a.(D$_)2014  0.994       0.995                 
## as.n((D$_)) -0.097      -0.105      -0.116     
## convergence code: 0
## boundary (singular) fit: see ?isSingular
anova(red_mod)
## Type III Analysis of Variance Table with Satterthwaite's method
##                                             Sum Sq Mean Sq NumDF DenDF F value
## redData$age_y                              0.53630 0.53630     1  1009 18.2001
## redData$feed_group                         0.03650 0.01825     2  1009  0.6193
## redData$genetic_group                      0.03650 0.03650     1  1009  1.2388
## as.factor(redData$sample_year)             0.55730 0.11146     5  1009  3.7825
## as.numeric(paste(redData$sample_interval)) 0.00115 0.00115     1  1009  0.0390
##                                               Pr(>F)    
## redData$age_y                              2.176e-05 ***
## redData$feed_group                          0.538542    
## redData$genetic_group                       0.265975    
## as.factor(redData$sample_year)              0.002125 ** 
## as.numeric(paste(redData$sample_interval))  0.843439    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
AIC(red_mod) 
## [1] -601.817

Reduced model of RLTL change with age as two level factor

red_mod_fac<-lmer(as.numeric(paste(redData$diff_rltl_res))~ redData$AgeGroup + redData$feed_group + redData$genetic_group+ as.factor(redData$sample_year) +  as.numeric(paste(redData$sample_interval))+(1|redData$recoded_id), na.action=na.exclude)
## boundary (singular) fit: see ?isSingular
summary(red_mod_fac)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: as.numeric(paste(redData$diff_rltl_res)) ~ redData$AgeGroup +  
##     redData$feed_group + redData$genetic_group + as.factor(redData$sample_year) +  
##     as.numeric(paste(redData$sample_interval)) + (1 | redData$recoded_id)
## 
## REML criterion at convergence: -627
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.6485 -0.6861  0.0272  0.6780  3.3620 
## 
## Random effects:
##  Groups             Name        Variance Std.Dev.
##  redData$recoded_id (Intercept) 0.00000  0.0000  
##  Residual                       0.02956  0.1719  
## Number of obs: 1020, groups:  redData$recoded_id, 305
## 
## Fixed effects:
##                                              Estimate Std. Error         df
## (Intercept)                                 2.988e-01  1.737e-01  1.009e+03
## redData$AgeGroupyoung                      -7.262e-02  1.872e-02  1.009e+03
## redData$feed_group2                        -5.832e-03  1.202e-02  1.009e+03
## redData$feed_groupNULL                     -1.786e-02  1.543e-02  1.009e+03
## redData$genetic_groupS                     -1.565e-02  1.131e-02  1.009e+03
## as.factor(redData$sample_year)2010         -3.503e-01  1.740e-01  1.009e+03
## as.factor(redData$sample_year)2011         -3.135e-01  1.737e-01  1.009e+03
## as.factor(redData$sample_year)2012         -2.886e-01  1.735e-01  1.009e+03
## as.factor(redData$sample_year)2013         -2.549e-01  1.737e-01  1.009e+03
## as.factor(redData$sample_year)2014         -2.564e-01  1.740e-01  1.009e+03
## as.numeric(paste(redData$sample_interval)) -7.172e-05  5.697e-05  1.009e+03
##                                            t value Pr(>|t|)    
## (Intercept)                                  1.720 0.085666 .  
## redData$AgeGroupyoung                       -3.880 0.000111 ***
## redData$feed_group2                         -0.485 0.627658    
## redData$feed_groupNULL                      -1.157 0.247431    
## redData$genetic_groupS                      -1.383 0.166882    
## as.factor(redData$sample_year)2010          -2.013 0.044340 *  
## as.factor(redData$sample_year)2011          -1.804 0.071475 .  
## as.factor(redData$sample_year)2012          -1.663 0.096544 .  
## as.factor(redData$sample_year)2013          -1.467 0.142608    
## as.factor(redData$sample_year)2014          -1.474 0.140820    
## as.numeric(paste(redData$sample_interval))  -1.259 0.208366    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) rdD$AG rdD$_2 rD$_NU rdD$_S a.(D$_)2010 a.(D$_)2011
## rdDt$AgGrpy -0.106                                                    
## rdDt$fd_gr2 -0.069  0.010                                             
## rdDt$f_NULL -0.034 -0.078  0.383                                      
## rdDt$gntc_S -0.060 -0.031 -0.020  0.237                               
## a.(D$_)2010 -0.983  0.009  0.029  0.013  0.024                        
## a.(D$_)2011 -0.987  0.033  0.030  0.005  0.022  0.991                 
## a.(D$_)2012 -0.989  0.044  0.031  0.010  0.024  0.992       0.995     
## a.(D$_)2013 -0.989  0.054  0.030  0.006  0.022  0.992       0.995     
## a.(D$_)2014 -0.987  0.051  0.027  0.006  0.019  0.991       0.995     
## as.n((D$_)) -0.050  0.404  0.044 -0.059  0.056 -0.088      -0.084     
##             a.(D$_)2012 a.(D$_)2013 a.(D$_)2014
## rdDt$AgGrpy                                    
## rdDt$fd_gr2                                    
## rdDt$f_NULL                                    
## rdDt$gntc_S                                    
## a.(D$_)2010                                    
## a.(D$_)2011                                    
## a.(D$_)2012                                    
## a.(D$_)2013  0.996                             
## a.(D$_)2014  0.996       0.996                 
## as.n((D$_)) -0.073      -0.078      -0.090     
## convergence code: 0
## boundary (singular) fit: see ?isSingular
anova(red_mod_fac)
## Type III Analysis of Variance Table with Satterthwaite's method
##                                             Sum Sq Mean Sq NumDF DenDF F value
## redData$AgeGroup                           0.44503 0.44503     1  1009 15.0564
## redData$feed_group                         0.03965 0.01982     2  1009  0.6707
## redData$genetic_group                      0.05656 0.05656     1  1009  1.9135
## as.factor(redData$sample_year)             0.89935 0.17987     5  1009  6.0855
## as.numeric(paste(redData$sample_interval)) 0.04684 0.04684     1  1009  1.5848
##                                               Pr(>F)    
## redData$AgeGroup                           0.0001111 ***
## redData$feed_group                         0.5115971    
## redData$genetic_group                      0.1668816    
## as.factor(redData$sample_year)              1.46e-05 ***
## as.numeric(paste(redData$sample_interval)) 0.2083656    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
AIC(red_mod_fac)
## [1] -600.9877
anova(red_mod, red_mod_fac)
## refitting model(s) with ML (instead of REML)
## Data: NULL
## Models:
## red_mod: as.numeric(paste(redData$diff_rltl_res)) ~ redData$age_y + redData$feed_group + 
## red_mod:     redData$genetic_group + as.factor(redData$sample_year) + 
## red_mod:     as.numeric(paste(redData$sample_interval)) + (1 | redData$recoded_id)
## red_mod_fac: as.numeric(paste(redData$diff_rltl_res)) ~ redData$AgeGroup + 
## red_mod_fac:     redData$feed_group + redData$genetic_group + as.factor(redData$sample_year) + 
## red_mod_fac:     as.numeric(paste(redData$sample_interval)) + (1 | redData$recoded_id)
##             npar     AIC     BIC logLik deviance Chisq Df Pr(>Chisq)
## red_mod       13 -685.40 -621.34 355.70  -711.40                    
## red_mod_fac   13 -682.27 -618.21 354.14  -708.27     0  0          1

Effect of milk productivity on RLTL change

Model of telomere change including milk productivity

#t test with average productivity in the model for animals that have productivity measurements:

PData<-subset(redData, redData$amp_305_lactation != "NULL")
nrow(PData)
## [1] 918
length(unique(PData$recoded_id))
## [1] 253
rescaled_prod<-as.numeric(paste(PData$amp_305_lactation))/1000
# milk productivity measures are on a completely differenc scale and are therefore re-scaled by dividing them by 1000

prod_mod<-lmer(as.numeric(paste(PData$diff_rltl_res))~ PData$age_y + PData$feed_group + PData$genetic_group+  as.factor(PData$sample_year) + as.numeric(paste(PData$sample_interval)) + rescaled_prod + (1|PData$recoded_id), na.action=na.exclude)
## boundary (singular) fit: see ?isSingular
summary(prod_mod)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: 
## as.numeric(paste(PData$diff_rltl_res)) ~ PData$age_y + PData$feed_group +  
##     PData$genetic_group + as.factor(PData$sample_year) + as.numeric(paste(PData$sample_interval)) +  
##     rescaled_prod + (1 | PData$recoded_id)
## 
## REML criterion at convergence: -555.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.5406 -0.6820  0.0102  0.7086  3.4856 
## 
## Random effects:
##  Groups           Name        Variance Std.Dev.
##  PData$recoded_id (Intercept) 0.00000  0.0000  
##  Residual                     0.02931  0.1712  
## Number of obs: 918, groups:  PData$recoded_id, 253
## 
## Fixed effects:
##                                            Estimate Std. Error         df
## (Intercept)                               1.954e-01  1.744e-01  9.060e+02
## PData$age_y                               2.069e-02  6.722e-03  9.060e+02
## PData$feed_group2                        -5.600e-03  1.305e-02  9.060e+02
## PData$feed_groupNULL                      9.214e-03  2.429e-02  9.060e+02
## PData$genetic_groupS                     -1.131e-02  1.222e-02  9.060e+02
## as.factor(PData$sample_year)2010         -3.381e-01  1.739e-01  9.060e+02
## as.factor(PData$sample_year)2011         -3.005e-01  1.737e-01  9.060e+02
## as.factor(PData$sample_year)2012         -2.783e-01  1.736e-01  9.060e+02
## as.factor(PData$sample_year)2013         -2.522e-01  1.745e-01  9.060e+02
## as.factor(PData$sample_year)2014         -2.746e-01  1.757e-01  9.060e+02
## as.numeric(paste(PData$sample_interval)) -2.573e-05  5.582e-05  9.060e+02
## rescaled_prod                             2.735e-03  2.840e-03  9.060e+02
##                                          t value Pr(>|t|)   
## (Intercept)                                1.121  0.26269   
## PData$age_y                                3.078  0.00215 **
## PData$feed_group2                         -0.429  0.66799   
## PData$feed_groupNULL                       0.379  0.70449   
## PData$genetic_groupS                      -0.926  0.35490   
## as.factor(PData$sample_year)2010          -1.945  0.05214 . 
## as.factor(PData$sample_year)2011          -1.730  0.08395 . 
## as.factor(PData$sample_year)2012          -1.603  0.10931   
## as.factor(PData$sample_year)2013          -1.446  0.14862   
## as.factor(PData$sample_year)2014          -1.563  0.11845   
## as.numeric(paste(PData$sample_interval))  -0.461  0.64498   
## rescaled_prod                              0.963  0.33571   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##              (Intr) PDt$g_ PDt$_2 PD$_NU PDt$_S a.(PD$_)2010 a.(PD$_)2011
## PData$age_y   0.033                                                      
## PDt$fd_grp2  -0.124 -0.163                                               
## PDt$fd_NULL  -0.127 -0.143  0.452                                        
## PDt$gntc_gS  -0.016  0.102 -0.130 -0.006                                 
## a.(PD$_)2010 -0.980 -0.049  0.050  0.044  0.010                          
## a.(PD$_)2011 -0.981 -0.069  0.054  0.044  0.006  0.991                   
## a.(PD$_)2012 -0.982 -0.087  0.057  0.051  0.005  0.991        0.994      
## a.(PD$_)2013 -0.979 -0.120  0.062  0.056  0.000  0.990        0.994      
## a.(PD$_)2014 -0.972 -0.161  0.065  0.060 -0.004  0.986        0.990      
## as.((PD$_))  -0.001  0.008  0.033 -0.064  0.059 -0.112       -0.115      
## rescald_prd  -0.164 -0.178  0.376  0.599 -0.283  0.055        0.057      
##              a.(PD$_)2012 a.(PD$_)2013 a.(PD$_)2014 a.((PD
## PData$age_y                                               
## PDt$fd_grp2                                               
## PDt$fd_NULL                                               
## PDt$gntc_gS                                               
## a.(PD$_)2010                                              
## a.(PD$_)2011                                              
## a.(PD$_)2012                                              
## a.(PD$_)2013  0.995                                       
## a.(PD$_)2014  0.993        0.995                          
## as.((PD$_))  -0.110       -0.121       -0.131             
## rescald_prd   0.060        0.068        0.071       -0.036
## convergence code: 0
## boundary (singular) fit: see ?isSingular
anova(prod_mod)
## Type III Analysis of Variance Table with Satterthwaite's method
##                                           Sum Sq  Mean Sq NumDF DenDF F value
## PData$age_y                              0.27757 0.277574     1   906  9.4712
## PData$feed_group                         0.01751 0.008753     2   906  0.2987
## PData$genetic_group                      0.02511 0.025108     1   906  0.8567
## as.factor(PData$sample_year)             0.49716 0.099433     5   906  3.3928
## as.numeric(paste(PData$sample_interval)) 0.00623 0.006226     1   906  0.2124
## rescaled_prod                            0.02719 0.027189     1   906  0.9277
##                                            Pr(>F)   
## PData$age_y                              0.002150 **
## PData$feed_group                         0.741877   
## PData$genetic_group                      0.354903   
## as.factor(PData$sample_year)             0.004824 **
## as.numeric(paste(PData$sample_interval)) 0.644976   
## rescaled_prod                            0.335712   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Analysis of change between first two measurements within first year and productive lifespan

Histogram of change in telomere length within the first year of life
ggplot(init_change_dat, aes(x=init_change_dat$ResDiff)) +geom_histogram(colour= "black", fill = mycoloursP[6])+ xlab("Initial RLTL change") + theme_classic(20)+  
  geom_vline(xintercept= mean(as.numeric(paste(init_change_dat$ResDiff))), colour="red") +
  annotate("text", x = -0.7, y = 28, label = "(*)", fontface =2, size= 10) +
    stat_function(
        fun = function(x, mean, sd, n){
             n/23 * dnorm(x = x, mean = mean, sd = sd)
        }, 
        args = with(init_change_dat, c(mean = mean(as.numeric(paste(init_change_dat$ResDiff))), 
                                 sd = sd(as.numeric(paste(init_change_dat$ResDiff))), 
                                 n= length(as.numeric(paste(init_change_dat$ResDiff))))))
## Warning: Use of `init_change_dat$ResDiff` is discouraged. Use `ResDiff` instead.

## Warning: Use of `init_change_dat$ResDiff` is discouraged. Use `ResDiff` instead.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Statistical normality tests
shapiro.test(init_change_dat$ResDiff)
## 
##  Shapiro-Wilk normality test
## 
## data:  init_change_dat$ResDiff
## W = 0.99401, p-value = 0.3057
ks.test(init_change_dat$ResDiff, pnorm)
## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  init_change_dat$ResDiff
## D = 0.39675, p-value < 2.2e-16
## alternative hypothesis: two-sided
qqnorm(init_change_dat$ResDiff)

Cox-proportional hazard model of productive lifespan
coxFit<-coxph(Surv(as.numeric(paste(init_change_dat$TimeMeasure)), init_change_dat$Event == 1)~init_change_dat$ResDiff)

summary(coxFit)
## Call:
## coxph(formula = Surv(as.numeric(paste(init_change_dat$TimeMeasure)), 
##     init_change_dat$Event == 1) ~ init_change_dat$ResDiff)
## 
##   n= 291, number of events= 227 
## 
##                            coef exp(coef) se(coef)      z Pr(>|z|)   
## init_change_dat$ResDiff -1.1406    0.3196   0.3914 -2.914  0.00356 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                         exp(coef) exp(-coef) lower .95 upper .95
## init_change_dat$ResDiff    0.3196      3.129    0.1484    0.6883
## 
## Concordance= 0.554  (se = 0.021 )
## Likelihood ratio test= 8.32  on 1 df,   p=0.004
## Wald test            = 8.49  on 1 df,   p=0.004
## Score (logrank) test = 8.48  on 1 df,   p=0.004
change telomere change measure to discrete scale to allow for plotting
ggplot(init_change_dat, aes(x=as.factor(init_change_dat$ResDiff_3G), y=init_change_dat$ResDiff)) + geom_boxplot() +xlab("RLTL change tertiles") + ylab("RLTL change")+ theme_classic(20)+ geom_jitter()
## Warning: Use of `init_change_dat$ResDiff_3G` is discouraged. Use `ResDiff_3G`
## instead.
## Warning: Use of `init_change_dat$ResDiff` is discouraged. Use `ResDiff` instead.
## Warning: Use of `init_change_dat$ResDiff_3G` is discouraged. Use `ResDiff_3G`
## instead.
## Warning: Use of `init_change_dat$ResDiff` is discouraged. Use `ResDiff` instead.

Cox proportional hazard model on discrete scale change measurements
coxFit<-coxph(Surv(init_change_dat$TimeMeasure, init_change_dat$Event == 1)~init_change_dat$ResDiff_3G)

summary(coxFit)
## Call:
## coxph(formula = Surv(init_change_dat$TimeMeasure, init_change_dat$Event == 
##     1) ~ init_change_dat$ResDiff_3G)
## 
##   n= 291, number of events= 227 
## 
##                               coef exp(coef) se(coef)      z Pr(>|z|)   
## init_change_dat$ResDiff_3G -0.2259    0.7978   0.0820 -2.755  0.00586 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                            exp(coef) exp(-coef) lower .95 upper .95
## init_change_dat$ResDiff_3G    0.7978      1.253    0.6793    0.9369
## 
## Concordance= 0.55  (se = 0.02 )
## Likelihood ratio test= 7.61  on 1 df,   p=0.006
## Wald test            = 7.59  on 1 df,   p=0.006
## Score (logrank) test = 7.65  on 1 df,   p=0.006
Plot discrete measures
#g1<- ggsurvplot(coxFit, data=init_change_dat, pval = 0.006, ggtheme = theme_minimal(), risk.table.y.text.col = T, risk.table.y.text = FALSE, legend.labs=c("Tertile 1 (most RLTL attrition)", "Tertile 2", "Tertile 3 (least RLTL attrition"), surv.median.line = "hv")+ guides(colour = guide_legend(nrow = 3))+ theme_survminer(40)

#g1

#ggsurvplot(coxFit, data=init_change_dat, pval = 0.006, ggtheme = theme_minimal(), risk.table.y.text.col = T, risk.table.y.text = FALSE, legend.labs=c("Tertile 1 (most RLTL attrition)", "Tertile 2", "Tertile 3 (least RLTL attrition"), surv.median.line = "hv")+theme_survminer(40) + xlab("Time in days")




#ggsave(file = "/Users/luise/Documents/Work/CattleTelomereLength/AgingCell/Images/InitialChange.pdf", print(g1))

Individual longitudinal RLTL data of animals with more than seven measurements

spaghettiData<-read.csv("/Users/luise/Documents/Work/CattleTelomereLength/LongTelomereChangePaper/git/TelomereChangeInDairyCattle/RC_SpaghettiData.csv")


ggplot(data = spaghettiData, aes(x = age_y, y = residual_rltl_t, group = recoded_id))+geom_point()+ geom_line()+facet_grid(.~recoded_id)+ 
  xlab("Age in years") + ylab("RLTL residuals")

ggplot(data = data, aes(x = sample_year, y = residual_rltl_t, group = recoded_id, colour= as.factor(data$recoded_id)))+geom_point()+ geom_line()+ 
  xlab("Age in years") + ylab("RLTL residuals") + 
  theme(legend.position = "none")
## Warning: Use of `data$recoded_id` is discouraged. Use `recoded_id` instead.

## Warning: Use of `data$recoded_id` is discouraged. Use `recoded_id` instead.

ggplot(data = data, aes(x = age_y , y = residual_rltl_t, group = recoded_id)) + geom_point(color= "white")+ geom_line(color= "white") + xlab("Age in years") + ylab("RLTL residuals") + theme(legend.position = "none") + geom_smooth(method = "lm", se = FALSE) + theme_classic()
## `geom_smooth()` using formula 'y ~ x'

# # average_resid rltl_change  abs_rltl_change MeanChangeRate MeanChangeRateAbs

deadD<-subset(data, data$herd_life != "NULL")
redDeadD<-subset(redData, redData$herd_life != "NULL")
uDeadData<-redDeadD[!duplicated(redDeadD$recoded_id),]

#1)
x<-15

AverRes<-uDeadData$average_resid
chF<-uDeadData$rltl_change

df1<-data.frame(AverRes, chF)


cor.test(AverRes, chF)
## 
##  Pearson's product-moment correlation
## 
## data:  AverRes and chF
## t = -3.0223, df = 239, p-value = 0.002782
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.31069811 -0.06712565
## sample estimates:
##        cor 
## -0.1918646
CorP1 <- ggplot(df1 , aes(x= AverRes, y=chF)) + geom_point(size= 2) +
  labs(x= "Mean RLTL", y= "Mean RLTL change") + stat_smooth(method=lm, formula = y~x)+
  theme_gray(x)
CorP1

#2) 
AchF<-uDeadData$abs_rltl_change

df2<-data.frame(AverRes, AchF)


cor.test(AverRes, AchF)
## 
##  Pearson's product-moment correlation
## 
## data:  AverRes and AchF
## t = 4.9383, df = 239, p-value = 1.48e-06
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.185034 0.414706
## sample estimates:
##       cor 
## 0.3042856
CorP2 <- ggplot(df2 , aes(x= AverRes, y=AchF)) + geom_point(size= 2) +
  labs(x= "Mean RLTL", y= "Mean absolute RLTL change") + stat_smooth(method=lm, formula = y~x)+
  theme_gray(x)
CorP2

#5) 

df5<-data.frame(chF, AchF)

cor.test(chF, AchF)
## 
##  Pearson's product-moment correlation
## 
## data:  chF and AchF
## t = -9.7941, df = 239, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.6196323 -0.4384552
## sample estimates:
##      cor 
## -0.53517
CorP5 <- ggplot(df5 , aes(x= chF, y=AchF)) + geom_point(size= 2) +
  labs(x= "Mean RLTL change", y= "Mean absolute RLTL change") +
  theme_gray(x)+ stat_smooth(method=lm, formula = y~x)+
  geom_vline(xintercept=0, linetype="dotted") +
  geom_abline(intercept=0, slope=-1, colour="red")+
  geom_abline(intercept=0, slope=1, colour= "blue")
CorP5

#grid.arrange(CorP1, CorP2, CorP3, CorP4, CorP5, CorP6, CorP7, CorP8, CorP9, CorP10, ncol=3 ) other plots can be found in 20170914_ChangeAnalysis

grid.arrange(CorP1, CorP2, CorP5, ncol=3 )

Correlation between subsquent measurements

cor.test(redData$residual_rltl_t,as.numeric(paste(redData$residual_rltl_tminus1)))
## 
##  Pearson's product-moment correlation
## 
## data:  redData$residual_rltl_t and as.numeric(paste(redData$residual_rltl_tminus1))
## t = 13.058, df = 1018, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.3249385 0.4301473
## sample estimates:
##       cor 
## 0.3787659
SubsPlot<-ggplot(redData, aes(x=as.numeric(paste(residual_rltl_tminus1)), y=redData$residual_rltl_t))+geom_point()+xlab("RLTL at time t-1")+ylab("RLTL at time t")+theme_gray(20)+geom_abline(intercept=0, slope=1, colour="red")+ stat_smooth(method=lm, formula=y~x)

SubsPlot
## Warning: Use of `redData$residual_rltl_t` is discouraged. Use `residual_rltl_t`
## instead.

## Warning: Use of `redData$residual_rltl_t` is discouraged. Use `residual_rltl_t`
## instead.

#Plot correlation between present and past change
redRedData<-subset(data, data$change_at_t_1!= "NULL")
nrow(redRedData)
## [1] 715
cor.test(as.numeric(paste(redRedData$diff_rltl_res)), as.numeric(paste(redRedData$change_at_t_1)))
## 
##  Pearson's product-moment correlation
## 
## data:  as.numeric(paste(redRedData$diff_rltl_res)) and as.numeric(paste(redRedData$change_at_t_1))
## t = -12.838, df = 713, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.4910345 -0.3718059
## sample estimates:
##        cor 
## -0.4333142
SubsPlot2<-ggplot(redRedData, aes(x=as.numeric(paste(redRedData$change_at_t_1)), y=as.numeric(paste(redRedData$diff_rltl_res))))+geom_point()+xlab("RLTL change at time t-1")+ylab("RLTL change at time t")+theme_gray(20)+geom_abline(intercept=0, slope=-1, colour="red")+ stat_smooth(method=lm, formula=y~x)

SubsPlot2
## Warning: Use of `redRedData$change_at_t_1` is discouraged. Use `change_at_t_1`
## instead.
## Warning: Use of `redRedData$diff_rltl_res` is discouraged. Use `diff_rltl_res`
## instead.
## Warning: Use of `redRedData$change_at_t_1` is discouraged. Use `change_at_t_1`
## instead.
## Warning: Use of `redRedData$diff_rltl_res` is discouraged. Use `diff_rltl_res`
## instead.

grid.arrange(SubsPlot, SubsPlot2, ncol=2)
## Warning: Use of `redData$residual_rltl_t` is discouraged. Use `residual_rltl_t`
## instead.
## Warning: Use of `redData$residual_rltl_t` is discouraged. Use `residual_rltl_t`
## instead.
## Warning: Use of `redRedData$change_at_t_1` is discouraged. Use `change_at_t_1`
## instead.
## Warning: Use of `redRedData$diff_rltl_res` is discouraged. Use `diff_rltl_res`
## instead.
## Warning: Use of `redRedData$change_at_t_1` is discouraged. Use `change_at_t_1`
## instead.
## Warning: Use of `redRedData$diff_rltl_res` is discouraged. Use `diff_rltl_res`
## instead.

#resRTLatTmin1<-as.numeric(paste(redData$residual_rltl_tminus1))

Productive lifespan

variable<- as.numeric(paste(uDeadData$herd_life))
#sqrt might improve normality test outcome slightly...

meanHL<-mean(as.numeric(paste(uDeadData$herd_life)))


prod_lsp_plot <- ggplot(uDeadData, aes(x = as.numeric(paste(uDeadData$herd_life)))) + 
    geom_histogram(colour= "black", fill = mycoloursP[4]) + xlab("Productive lifespan in days") +  
  geom_vline(xintercept= mean(as.numeric(paste(uDeadData$herd_life))), colour="red")+
  annotate("text", x = 0, y = 17.5, label = "(D)", fontface =2, size= 10) + theme_classic(20) +
    stat_function(
        fun = function(x, mean, sd, n){
             n*120 * dnorm(x = x, mean = mean, sd = sd)
        }, 
        args = with(uDeadData, c(mean = mean(as.numeric(paste(uDeadData$herd_life))), 
                                 sd = sd(as.numeric(paste(uDeadData$herd_life))), 
                                 n= length(as.numeric(paste(uDeadData$herd_life))))))


prod_lsp_plot
## Warning: Use of `uDeadData$herd_life` is discouraged. Use `herd_life` instead.

## Warning: Use of `uDeadData$herd_life` is discouraged. Use `herd_life` instead.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

# statistical tests for normal distribution of RTL measurements
shapiro.test(variable)
## 
##  Shapiro-Wilk normality test
## 
## data:  variable
## W = 0.98408, p-value = 0.008487
ks.test(variable, pnorm)
## Warning in ks.test(variable, pnorm): ties should not be present for the
## Kolmogorov-Smirnov test
## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  variable
## D = 1, p-value < 2.2e-16
## alternative hypothesis: two-sided
qqnorm(variable)

Assemble plots for supplemenrary file

grid.arrange(hist1, hist2, hist3, prod_lsp_plot, ncol=2)
## Warning: Use of `uData$SamplesPerAnimal` is discouraged. Use `SamplesPerAnimal`
## instead.
## Warning: Use of `data$residual_rltl_t` is discouraged. Use `residual_rltl_t`
## instead.

## Warning: Use of `data$residual_rltl_t` is discouraged. Use `residual_rltl_t`
## instead.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning: Use of `redData$diff_rltl_res` is discouraged. Use `diff_rltl_res`
## instead.
## Warning: Use of `redData$diff_rltl_res` is discouraged. Use `diff_rltl_res`
## instead.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning: Use of `uDeadData$herd_life` is discouraged. Use `herd_life` instead.
## Warning: Use of `uDeadData$herd_life` is discouraged. Use `herd_life` instead.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Histograms of (A) Number of RLTL measurements per animal, (B) RLTL residuals after correcting for qPCR plate and row, (C) measurement to measurement change in RLTL residuals (after correcting for qPCR plate and row) and (D) productive lifespan in days; Red lines indicate mean.

Histograms of (A) Number of RLTL measurements per animal, (B) RLTL residuals after correcting for qPCR plate and row, (C) measurement to measurement change in RLTL residuals (after correcting for qPCR plate and row) and (D) productive lifespan in days; Red lines indicate mean.

Backwards elimination of all non-significant effects

#Age in years as factor with reference year being 0 (Age at t-1)
mod1f<-lmer(as.numeric(paste(redData$diff_rltl_res))~ as.factor(redData$age_y) + (1|redData$recoded_id), na.action=na.exclude)
## boundary (singular) fit: see ?isSingular
summary(mod1f)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: as.numeric(paste(redData$diff_rltl_res)) ~ as.factor(redData$age_y) +  
##     (1 | redData$recoded_id)
## 
## REML criterion at convergence: -691.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.6127 -0.6744  0.0379  0.6816  3.5134 
## 
## Random effects:
##  Groups             Name        Variance Std.Dev.
##  redData$recoded_id (Intercept) 0.00000  0.0000  
##  Residual                       0.02871  0.1694  
## Number of obs: 1020, groups:  redData$recoded_id, 305
## 
## Fixed effects:
##                             Estimate Std. Error         df t value Pr(>|t|)  
## (Intercept)                 -0.10057    0.06404 1013.00000  -1.570   0.1167  
## as.factor(redData$age_y)1   -0.01458    0.06482 1013.00000  -0.225   0.8221  
## as.factor(redData$age_y)2    0.09171    0.06484 1013.00000   1.414   0.1575  
## as.factor(redData$age_y)3    0.11267    0.06503 1013.00000   1.733   0.0835 .
## as.factor(redData$age_y)4    0.09686    0.06558 1013.00000   1.477   0.1400  
## as.factor(redData$age_y)5    0.12766    0.06800 1013.00000   1.877   0.0608 .
## as.factor(redData$age_y)6    0.07425    0.07441 1013.00000   0.998   0.3186  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) a.(D$_)1 a.(D$_)2 a.(D$_)3 a.(D$_)4 a.(D$_)5
## as.fc(D$_)1 -0.988                                             
## as.fc(D$_)2 -0.988  0.976                                      
## as.fc(D$_)3 -0.985  0.973    0.973                             
## as.fc(D$_)4 -0.977  0.965    0.965    0.962                    
## as.fc(D$_)5 -0.942  0.931    0.930    0.928    0.920           
## as.fc(D$_)6 -0.861  0.850    0.850    0.848    0.840    0.811  
## convergence code: 0
## boundary (singular) fit: see ?isSingular
anova(mod1f)
## Type III Analysis of Variance Table with Satterthwaite's method
##                          Sum Sq Mean Sq NumDF DenDF F value    Pr(>F)    
## as.factor(redData$age_y) 2.9149 0.48582     6  1013  16.921 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
################
################
#Test

#mod1f<-lmer(abs(as.numeric(paste(redData$diff_rltl_res)))~ 1 + (1|redData$recoded_id), na.action=na.exclude)

#summary(mod1f)
#anova(mod1f)


0.0006458/(0.0006458+0.0111165)
## [1] 0.05490423
##############
###############

redData$R_age_y<-as.factor(redData$age_y)

#lm without random effect to do likelihood ratio test

modLM<-glm(as.numeric(paste(redData$diff_rltl_res))~ as.factor(redData$age_y))
anova(mod1f, modLM)
## refitting model(s) with ML (instead of REML)
## Data: NULL
## Models:
## modLM: as.numeric(paste(redData$diff_rltl_res)) ~ as.factor(redData$age_y)
## mod1f: as.numeric(paste(redData$diff_rltl_res)) ~ as.factor(redData$age_y) + 
## mod1f:     (1 | redData$recoded_id)
##       npar     AIC     BIC logLik deviance Chisq Df Pr(>Chisq)
## modLM    8 -717.85 -678.43 366.93  -733.85                    
## mod1f    9 -715.85 -671.50 366.93  -733.85     0  1          1
logLik(mod1f)
## 'log Lik.' 345.6007 (df=9)
logLik(modLM)
## 'log Lik.' 366.9254 (df=8)
#Age in years as factor with reference year being 1 (Age at t-1)
contrasts(redData$R_age_y) <- contr.treatment(levels(redData$R_age_y), base=which(levels(redData$R_age_y) == '1'))

mod1f<-lmer(as.numeric(paste(redData$diff_rltl_res))~ as.factor(redData$R_age_y)+ as.factor(redData$sample_year) + as.numeric(paste(redData$sample_interval)) + (1|redData$recoded_id), na.action=na.exclude)
## boundary (singular) fit: see ?isSingular
summary(mod1f)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: 
## as.numeric(paste(redData$diff_rltl_res)) ~ as.factor(redData$R_age_y) +  
##     as.factor(redData$sample_year) + as.numeric(paste(redData$sample_interval)) +  
##     (1 | redData$recoded_id)
## 
## REML criterion at convergence: -653.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.5347 -0.6741  0.0278  0.6901  3.5959 
## 
## Random effects:
##  Groups             Name        Variance Std.Dev.
##  redData$recoded_id (Intercept) 0.00000  0.0000  
##  Residual                       0.02868  0.1693  
## Number of obs: 1020, groups:  redData$recoded_id, 305
## 
## Fixed effects:
##                                              Estimate Std. Error         df
## (Intercept)                                 2.363e-01  1.842e-01  1.007e+03
## as.factor(redData$R_age_y)0                -3.155e-02  7.258e-02  1.007e+03
## as.factor(redData$R_age_y)2                 9.681e-02  1.561e-02  1.007e+03
## as.factor(redData$R_age_y)3                 1.131e-01  1.827e-02  1.007e+03
## as.factor(redData$R_age_y)4                 9.575e-02  2.220e-02  1.007e+03
## as.factor(redData$R_age_y)5                 1.322e-01  3.037e-02  1.007e+03
## as.factor(redData$R_age_y)6                 7.947e-02  4.352e-02  1.007e+03
## as.factor(redData$sample_year)2010         -3.634e-01  1.831e-01  1.007e+03
## as.factor(redData$sample_year)2011         -3.539e-01  1.839e-01  1.007e+03
## as.factor(redData$sample_year)2012         -3.420e-01  1.838e-01  1.007e+03
## as.factor(redData$sample_year)2013         -3.257e-01  1.838e-01  1.007e+03
## as.factor(redData$sample_year)2014         -3.409e-01  1.839e-01  1.007e+03
## as.numeric(paste(redData$sample_interval)) -3.703e-06  5.319e-05  1.007e+03
##                                            t value Pr(>|t|)    
## (Intercept)                                  1.282   0.2000    
## as.factor(redData$R_age_y)0                 -0.435   0.6639    
## as.factor(redData$R_age_y)2                  6.201 8.19e-10 ***
## as.factor(redData$R_age_y)3                  6.192 8.64e-10 ***
## as.factor(redData$R_age_y)4                  4.314 1.76e-05 ***
## as.factor(redData$R_age_y)5                  4.352 1.49e-05 ***
## as.factor(redData$R_age_y)6                  1.826   0.0681 .  
## as.factor(redData$sample_year)2010          -1.984   0.0475 *  
## as.factor(redData$sample_year)2011          -1.924   0.0546 .  
## as.factor(redData$sample_year)2012          -1.861   0.0631 .  
## as.factor(redData$sample_year)2013          -1.772   0.0768 .  
## as.factor(redData$sample_year)2014          -1.854   0.0641 .  
## as.numeric(paste(redData$sample_interval))  -0.070   0.9445    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation matrix not shown by default, as p = 13 > 12.
## Use print(x, correlation=TRUE)  or
##     vcov(x)        if you need it
## convergence code: 0
## boundary (singular) fit: see ?isSingular
anova(mod1f)
## Type III Analysis of Variance Table with Satterthwaite's method
##                                             Sum Sq  Mean Sq NumDF DenDF F value
## as.factor(redData$R_age_y)                 1.49717 0.249528     6  1007  8.7016
## as.factor(redData$sample_year)             0.20720 0.041441     5  1007  1.4451
## as.numeric(paste(redData$sample_interval)) 0.00014 0.000139     1  1007  0.0048
##                                               Pr(>F)    
## as.factor(redData$R_age_y)                 2.959e-09 ***
## as.factor(redData$sample_year)                0.2054    
## as.numeric(paste(redData$sample_interval))    0.9445    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
exclude2009<-subset(redData, redData$sample_year != "2009")

mod1f<-lmer(as.numeric(paste(exclude2009$diff_rltl_res))~ as.factor(exclude2009$R_age_y)+ as.factor(exclude2009$sample_year) + as.numeric(paste(exclude2009$sample_interval)) + (1|exclude2009$recoded_id), na.action=na.exclude)
## boundary (singular) fit: see ?isSingular
summary(mod1f)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: 
## as.numeric(paste(exclude2009$diff_rltl_res)) ~ as.factor(exclude2009$R_age_y) +  
##     as.factor(exclude2009$sample_year) + as.numeric(paste(exclude2009$sample_interval)) +  
##     (1 | exclude2009$recoded_id)
## 
## REML criterion at convergence: -653.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.5347 -0.6744  0.0286  0.6916  3.5959 
## 
## Random effects:
##  Groups                 Name        Variance Std.Dev.
##  exclude2009$recoded_id (Intercept) 0.00000  0.0000  
##  Residual                           0.02868  0.1693  
## Number of obs: 1019, groups:  exclude2009$recoded_id, 305
## 
## Fixed effects:
##                                                  Estimate Std. Error         df
## (Intercept)                                    -1.272e-01  2.672e-02  1.007e+03
## as.factor(exclude2009$R_age_y)0                -3.155e-02  7.258e-02  1.007e+03
## as.factor(exclude2009$R_age_y)2                 9.681e-02  1.561e-02  1.007e+03
## as.factor(exclude2009$R_age_y)3                 1.131e-01  1.827e-02  1.007e+03
## as.factor(exclude2009$R_age_y)4                 9.575e-02  2.220e-02  1.007e+03
## as.factor(exclude2009$R_age_y)5                 1.322e-01  3.037e-02  1.007e+03
## as.factor(exclude2009$R_age_y)6                 7.947e-02  4.352e-02  1.007e+03
## as.factor(exclude2009$sample_year)2011          9.552e-03  2.315e-02  1.007e+03
## as.factor(exclude2009$sample_year)2012          2.146e-02  2.309e-02  1.007e+03
## as.factor(exclude2009$sample_year)2013          3.772e-02  2.459e-02  1.007e+03
## as.factor(exclude2009$sample_year)2014          2.254e-02  2.724e-02  1.007e+03
## as.numeric(paste(exclude2009$sample_interval)) -3.703e-06  5.319e-05  1.007e+03
##                                                t value Pr(>|t|)    
## (Intercept)                                     -4.759 2.23e-06 ***
## as.factor(exclude2009$R_age_y)0                 -0.435   0.6639    
## as.factor(exclude2009$R_age_y)2                  6.201 8.19e-10 ***
## as.factor(exclude2009$R_age_y)3                  6.192 8.64e-10 ***
## as.factor(exclude2009$R_age_y)4                  4.314 1.76e-05 ***
## as.factor(exclude2009$R_age_y)5                  4.352 1.49e-05 ***
## as.factor(exclude2009$R_age_y)6                  1.826   0.0681 .  
## as.factor(exclude2009$sample_year)2011           0.413   0.6799    
## as.factor(exclude2009$sample_year)2012           0.930   0.3528    
## as.factor(exclude2009$sample_year)2013           1.534   0.1253    
## as.factor(exclude2009$sample_year)2014           0.827   0.4083    
## as.numeric(paste(exclude2009$sample_interval))  -0.070   0.9445    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##                (Intr) a.(2009$R__)0 a.(2009$R__)2 a.(2009$R__)3 a.(2009$R__)4
## a.(2009$R__)0  -0.287                                                        
## a.(2009$R__)2  -0.031  0.057                                                 
## a.(2009$R__)3  -0.007  0.064         0.560                                   
## a.(2009$R__)4   0.019  0.060         0.495         0.585                     
## a.(2009$R__)5  -0.015  0.066         0.362         0.456         0.483       
## a.(2009$R__)6  -0.107  0.078         0.251         0.318         0.335       
## a.(2009$_)2011 -0.520  0.113        -0.267        -0.162        -0.133       
## a.(2009$_)2012 -0.566  0.117        -0.329        -0.336        -0.236       
## a.(2009$_)2013 -0.486  0.081        -0.405        -0.453        -0.433       
## a.(2009$_)2014 -0.362  0.013        -0.362        -0.484        -0.516       
## a.((2009$_)    -0.726  0.232        -0.014        -0.028        -0.062       
##                a.(2009$R__)5 a.(2009$R__)6 a.(2009$_)2011 a.(2009$_)2012
## a.(2009$R__)0                                                           
## a.(2009$R__)2                                                           
## a.(2009$R__)3                                                           
## a.(2009$R__)4                                                           
## a.(2009$R__)5                                                           
## a.(2009$R__)6   0.313                                                   
## a.(2009$_)2011 -0.097        -0.070                                     
## a.(2009$_)2012 -0.173        -0.115         0.719                       
## a.(2009$_)2013 -0.281        -0.192         0.702          0.775        
## a.(2009$_)2014 -0.480        -0.362         0.632          0.708        
## a.((2009$_)    -0.008         0.127        -0.017          0.049        
##                a.(2009$_)2013 a.(2009$_)2014
## a.(2009$R__)0                               
## a.(2009$R__)2                               
## a.(2009$R__)3                               
## a.(2009$R__)4                               
## a.(2009$R__)5                               
## a.(2009$R__)6                               
## a.(2009$_)2011                              
## a.(2009$_)2012                              
## a.(2009$_)2013                              
## a.(2009$_)2014  0.790                       
## a.((2009$_)    -0.010         -0.108        
## convergence code: 0
## boundary (singular) fit: see ?isSingular
anova(mod1f)
## Type III Analysis of Variance Table with Satterthwaite's method
##                                                 Sum Sq  Mean Sq NumDF DenDF
## as.factor(exclude2009$R_age_y)                 1.49717 0.249528     6  1007
## as.factor(exclude2009$sample_year)             0.09885 0.024712     4  1007
## as.numeric(paste(exclude2009$sample_interval)) 0.00014 0.000139     1  1007
##                                                F value    Pr(>F)    
## as.factor(exclude2009$R_age_y)                  8.7016 2.959e-09 ***
## as.factor(exclude2009$sample_year)              0.8618    0.4864    
## as.numeric(paste(exclude2009$sample_interval))  0.0048    0.9445    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
df<-as.data.frame(round(coef(summary(mod1f)),3))
#write.csv(df, "/Users/luiseseeker/Desktop/PNAS/S2_Table_AgeFactor.csv")

#lsmeans(mod1f,pairwise~as.factor(redData$R_age_y), adjust="Bonferroni") 


#Age in years as factor with reference year being 2 (Age at t-1)
contrasts(redData$R_age_y) <- contr.treatment(levels(redData$R_age_y), base=which(levels(redData$R_age_y) == '2'))

mod1f<-lmer(as.numeric(paste(redData$diff_rltl_res))~ as.factor(redData$R_age_y) + (1|redData$recoded_id), na.action=na.exclude)
## boundary (singular) fit: see ?isSingular
summary(mod1f)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: 
## as.numeric(paste(redData$diff_rltl_res)) ~ as.factor(redData$R_age_y) +  
##     (1 | redData$recoded_id)
## 
## REML criterion at convergence: -691.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.6127 -0.6744  0.0379  0.6816  3.5134 
## 
## Random effects:
##  Groups             Name        Variance Std.Dev.
##  redData$recoded_id (Intercept) 0.00000  0.0000  
##  Residual                       0.02871  0.1694  
## Number of obs: 1020, groups:  redData$recoded_id, 305
## 
## Fixed effects:
##                               Estimate Std. Error         df t value Pr(>|t|)
## (Intercept)                 -8.857e-03  1.011e-02  1.013e+03  -0.876    0.381
## as.factor(redData$R_age_y)0 -9.171e-02  6.484e-02  1.013e+03  -1.414    0.158
## as.factor(redData$R_age_y)1 -1.063e-01  1.422e-02  1.013e+03  -7.474 1.67e-13
## as.factor(redData$R_age_y)3  2.096e-02  1.514e-02  1.013e+03   1.385    0.166
## as.factor(redData$R_age_y)4  5.150e-03  1.737e-02  1.013e+03   0.297    0.767
## as.factor(redData$R_age_y)5  3.595e-02  2.498e-02  1.013e+03   1.439    0.150
## as.factor(redData$R_age_y)6 -1.746e-02  3.921e-02  1.013e+03  -0.445    0.656
##                                
## (Intercept)                    
## as.factor(redData$R_age_y)0    
## as.factor(redData$R_age_y)1 ***
## as.factor(redData$R_age_y)3    
## as.factor(redData$R_age_y)4    
## as.factor(redData$R_age_y)5    
## as.factor(redData$R_age_y)6    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) a.(D$R__)0 a.(D$R__)1 a.(D$R__)3 a.(D$R__)4 a.(D$R__)5
## as.(D$R__)0 -0.156                                                       
## as.(D$R__)1 -0.711  0.111                                                
## as.(D$R__)3 -0.668  0.104      0.475                                     
## as.(D$R__)4 -0.582  0.091      0.414      0.389                          
## as.(D$R__)5 -0.405  0.063      0.288      0.270      0.236               
## as.(D$R__)6 -0.258  0.040      0.183      0.172      0.150      0.104    
## convergence code: 0
## boundary (singular) fit: see ?isSingular
anova(mod1f)
## Type III Analysis of Variance Table with Satterthwaite's method
##                            Sum Sq Mean Sq NumDF DenDF F value    Pr(>F)    
## as.factor(redData$R_age_y) 2.9149 0.48582     6  1013  16.921 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Age in years as factor with reference year being 3 (Age at t-1)
contrasts(redData$R_age_y) <- contr.treatment(levels(redData$R_age_y), base=which(levels(redData$R_age_y) == '3'))

mod1f<-lmer(as.numeric(paste(redData$diff_rltl_res))~ as.factor(redData$R_age_y) + (1|redData$recoded_id), na.action=na.exclude)
## boundary (singular) fit: see ?isSingular
summary(mod1f)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: 
## as.numeric(paste(redData$diff_rltl_res)) ~ as.factor(redData$R_age_y) +  
##     (1 | redData$recoded_id)
## 
## REML criterion at convergence: -691.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.6127 -0.6744  0.0379  0.6816  3.5134 
## 
## Random effects:
##  Groups             Name        Variance Std.Dev.
##  redData$recoded_id (Intercept) 0.00000  0.0000  
##  Residual                       0.02871  0.1694  
## Number of obs: 1020, groups:  redData$recoded_id, 305
## 
## Fixed effects:
##                               Estimate Std. Error         df t value Pr(>|t|)
## (Intercept)                    0.01211    0.01127 1013.00000   1.074   0.2830
## as.factor(redData$R_age_y)0   -0.11267    0.06503 1013.00000  -1.733   0.0835
## as.factor(redData$R_age_y)1   -0.12725    0.01507 1013.00000  -8.444   <2e-16
## as.factor(redData$R_age_y)2   -0.02096    0.01514 1013.00000  -1.385   0.1665
## as.factor(redData$R_age_y)4   -0.01581    0.01807 1013.00000  -0.875   0.3816
## as.factor(redData$R_age_y)5    0.01499    0.02548 1013.00000   0.588   0.5565
## as.factor(redData$R_age_y)6   -0.03843    0.03953 1013.00000  -0.972   0.3312
##                                
## (Intercept)                    
## as.factor(redData$R_age_y)0 .  
## as.factor(redData$R_age_y)1 ***
## as.factor(redData$R_age_y)2    
## as.factor(redData$R_age_y)4    
## as.factor(redData$R_age_y)5    
## as.factor(redData$R_age_y)6    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) a.(D$R__)0 a.(D$R__)1 a.(D$R__)2 a.(D$R__)4 a.(D$R__)5
## as.(D$R__)0 -0.173                                                       
## as.(D$R__)1 -0.748  0.130                                                
## as.(D$R__)2 -0.744  0.129      0.557                                     
## as.(D$R__)4 -0.624  0.108      0.467      0.464                          
## as.(D$R__)5 -0.442  0.077      0.331      0.329      0.276               
## as.(D$R__)6 -0.285  0.049      0.213      0.212      0.178      0.126    
## convergence code: 0
## boundary (singular) fit: see ?isSingular
anova(mod1f)
## Type III Analysis of Variance Table with Satterthwaite's method
##                            Sum Sq Mean Sq NumDF DenDF F value    Pr(>F)    
## as.factor(redData$R_age_y) 2.9149 0.48582     6  1013  16.921 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Age in years as factor with reference year being 4 (Age at t-1)
contrasts(redData$R_age_y) <- contr.treatment(levels(redData$R_age_y), base=which(levels(redData$R_age_y) == '4'))

mod1f<-lmer(as.numeric(paste(redData$diff_rltl_res))~ as.factor(redData$R_age_y) + (1|redData$recoded_id), na.action=na.exclude)
## boundary (singular) fit: see ?isSingular
summary(mod1f)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: 
## as.numeric(paste(redData$diff_rltl_res)) ~ as.factor(redData$R_age_y) +  
##     (1 | redData$recoded_id)
## 
## REML criterion at convergence: -691.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.6127 -0.6744  0.0379  0.6816  3.5134 
## 
## Random effects:
##  Groups             Name        Variance Std.Dev.
##  redData$recoded_id (Intercept) 0.00000  0.0000  
##  Residual                       0.02871  0.1694  
## Number of obs: 1020, groups:  redData$recoded_id, 305
## 
## Fixed effects:
##                               Estimate Std. Error         df t value Pr(>|t|)
## (Intercept)                 -3.707e-03  1.412e-02  1.013e+03  -0.263    0.793
## as.factor(redData$R_age_y)0 -9.686e-02  6.558e-02  1.013e+03  -1.477    0.140
## as.factor(redData$R_age_y)1 -1.114e-01  1.730e-02  1.013e+03  -6.440 1.84e-10
## as.factor(redData$R_age_y)2 -5.150e-03  1.737e-02  1.013e+03  -0.297    0.767
## as.factor(redData$R_age_y)3  1.581e-02  1.807e-02  1.013e+03   0.875    0.382
## as.factor(redData$R_age_y)5  3.080e-02  2.686e-02  1.013e+03   1.147    0.252
## as.factor(redData$R_age_y)6 -2.261e-02  4.043e-02  1.013e+03  -0.559    0.576
##                                
## (Intercept)                    
## as.factor(redData$R_age_y)0    
## as.factor(redData$R_age_y)1 ***
## as.factor(redData$R_age_y)2    
## as.factor(redData$R_age_y)3    
## as.factor(redData$R_age_y)5    
## as.factor(redData$R_age_y)6    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) a.(D$R__)0 a.(D$R__)1 a.(D$R__)2 a.(D$R__)3 a.(D$R__)5
## as.(D$R__)0 -0.215                                                       
## as.(D$R__)1 -0.816  0.176                                                
## as.(D$R__)2 -0.813  0.175      0.664                                     
## as.(D$R__)3 -0.782  0.168      0.638      0.635                          
## as.(D$R__)5 -0.526  0.113      0.429      0.427      0.411               
## as.(D$R__)6 -0.349  0.075      0.285      0.284      0.273      0.184    
## convergence code: 0
## boundary (singular) fit: see ?isSingular
anova(mod1f)
## Type III Analysis of Variance Table with Satterthwaite's method
##                            Sum Sq Mean Sq NumDF DenDF F value    Pr(>F)    
## as.factor(redData$R_age_y) 2.9149 0.48582     6  1013  16.921 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Age in years as factor with reference year being 5 (Age at t-1)
contrasts(redData$R_age_y) <- contr.treatment(levels(redData$R_age_y), base=which(levels(redData$R_age_y) == '5'))

mod1f<-lmer(as.numeric(paste(redData$diff_rltl_res))~ as.factor(redData$R_age_y) + (1|redData$recoded_id), na.action=na.exclude)
## boundary (singular) fit: see ?isSingular
summary(mod1f)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: 
## as.numeric(paste(redData$diff_rltl_res)) ~ as.factor(redData$R_age_y) +  
##     (1 | redData$recoded_id)
## 
## REML criterion at convergence: -691.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.6127 -0.6744  0.0379  0.6816  3.5134 
## 
## Random effects:
##  Groups             Name        Variance Std.Dev.
##  redData$recoded_id (Intercept) 0.00000  0.0000  
##  Residual                       0.02871  0.1694  
## Number of obs: 1020, groups:  redData$recoded_id, 305
## 
## Fixed effects:
##                               Estimate Std. Error         df t value Pr(>|t|)
## (Intercept)                    0.02709    0.02285 1013.00000   1.186   0.2360
## as.factor(redData$R_age_y)0   -0.12766    0.06800 1013.00000  -1.877   0.0608
## as.factor(redData$R_age_y)1   -0.14224    0.02494 1013.00000  -5.703 1.55e-08
## as.factor(redData$R_age_y)2   -0.03595    0.02498 1013.00000  -1.439   0.1505
## as.factor(redData$R_age_y)3   -0.01499    0.02548 1013.00000  -0.588   0.5565
## as.factor(redData$R_age_y)4   -0.03080    0.02686 1013.00000  -1.147   0.2518
## as.factor(redData$R_age_y)6   -0.05341    0.04424 1013.00000  -1.207   0.2276
##                                
## (Intercept)                    
## as.factor(redData$R_age_y)0 .  
## as.factor(redData$R_age_y)1 ***
## as.factor(redData$R_age_y)2    
## as.factor(redData$R_age_y)3    
## as.factor(redData$R_age_y)4    
## as.factor(redData$R_age_y)6    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) a.(D$R__)0 a.(D$R__)1 a.(D$R__)2 a.(D$R__)3 a.(D$R__)4
## as.(D$R__)0 -0.336                                                       
## as.(D$R__)1 -0.916  0.308                                                
## as.(D$R__)2 -0.914  0.307      0.838                                     
## as.(D$R__)3 -0.897  0.301      0.822      0.820                          
## as.(D$R__)4 -0.851  0.286      0.779      0.778      0.763               
## as.(D$R__)6 -0.516  0.174      0.473      0.472      0.463      0.439    
## convergence code: 0
## boundary (singular) fit: see ?isSingular
anova(mod1f)
## Type III Analysis of Variance Table with Satterthwaite's method
##                            Sum Sq Mean Sq NumDF DenDF F value    Pr(>F)    
## as.factor(redData$R_age_y) 2.9149 0.48582     6  1013  16.921 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
lm<-lm(as.numeric(paste(redData$diff_rltl_res))~as.factor(redData$age_y))
summary(lm)
## 
## Call:
## lm(formula = as.numeric(paste(redData$diff_rltl_res)) ~ as.factor(redData$age_y))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.61215 -0.11427  0.00643  0.11550  0.59533 
## 
## Coefficients:
##                           Estimate Std. Error t value Pr(>|t|)  
## (Intercept)               -0.10057    0.06404  -1.570   0.1167  
## as.factor(redData$age_y)1 -0.01458    0.06482  -0.225   0.8221  
## as.factor(redData$age_y)2  0.09171    0.06484   1.414   0.1575  
## as.factor(redData$age_y)3  0.11267    0.06503   1.733   0.0835 .
## as.factor(redData$age_y)4  0.09686    0.06558   1.477   0.1400  
## as.factor(redData$age_y)5  0.12766    0.06800   1.877   0.0608 .
## as.factor(redData$age_y)6  0.07425    0.07441   0.998   0.3186  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1694 on 1013 degrees of freedom
## Multiple R-squared:  0.09109,    Adjusted R-squared:  0.08571 
## F-statistic: 16.92 on 6 and 1013 DF,  p-value: < 2.2e-16
anova(lm)
## Analysis of Variance Table
## 
## Response: as.numeric(paste(redData$diff_rltl_res))
##                            Df  Sum Sq Mean Sq F value    Pr(>F)    
## as.factor(redData$age_y)    6  2.9149 0.48582  16.921 < 2.2e-16 ***
## Residuals                1013 29.0849 0.02871                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nrow(subset(redData, redData$age_y==6))
## [1] 20
#lsmeans(lm,pairwise~as.factor(redData$age_y), adjust="Bonferroni")

Test if RLTL depends on RLTL measurement at t-1

mod1<-lmer(as.numeric(paste(redData$residual_rltl_t))~ as.numeric(paste(redData$residual_rltl_tminus1))+redData$age_y + redData$feed_group + redData$genetic_group+ redData$birth_year+ as.numeric(paste(redData$sample_interval))+ redData$health_event_2week+ (1|redData$recoded_id), na.action=na.exclude)

summary(mod1)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: 
## as.numeric(paste(redData$residual_rltl_t)) ~ as.numeric(paste(redData$residual_rltl_tminus1)) +  
##     redData$age_y + redData$feed_group + redData$genetic_group +  
##     redData$birth_year + as.numeric(paste(redData$sample_interval)) +  
##     redData$health_event_2week + (1 | redData$recoded_id)
## 
## REML criterion at convergence: -1064.1
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.51548 -0.66303 -0.06712  0.55208  3.08068 
## 
## Random effects:
##  Groups             Name        Variance Std.Dev.
##  redData$recoded_id (Intercept) 0.005512 0.07425 
##  Residual                       0.015346 0.12388 
## Number of obs: 1020, groups:  redData$recoded_id, 305
## 
## Fixed effects:
##                                                    Estimate Std. Error
## (Intercept)                                       1.588e+01  1.126e+01
## as.numeric(paste(redData$residual_rltl_tminus1))  1.577e-01  2.921e-02
## redData$age_y                                     6.145e-03  3.474e-03
## redData$feed_group2                               6.969e-03  1.336e-02
## redData$feed_groupNULL                            8.637e-03  1.640e-02
## redData$genetic_groupS                            9.229e-03  1.239e-02
## redData$birth_year                               -7.928e-03  5.603e-03
## as.numeric(paste(redData$sample_interval))        1.349e-05  4.020e-05
## redData$health_event_2week                       -4.652e-03  1.920e-02
##                                                          df t value Pr(>|t|)
## (Intercept)                                       2.079e+02   1.410   0.1601
## as.numeric(paste(redData$residual_rltl_tminus1))  9.784e+02   5.397 8.49e-08
## redData$age_y                                     6.804e+02   1.769   0.0773
## redData$feed_group2                               1.446e+02   0.522   0.6027
## redData$feed_groupNULL                            1.787e+02   0.527   0.5990
## redData$genetic_groupS                            1.531e+02   0.745   0.4573
## redData$birth_year                                2.078e+02  -1.415   0.1586
## as.numeric(paste(redData$sample_interval))        8.790e+02   0.335   0.7373
## redData$health_event_2week                        8.576e+02  -0.242   0.8086
##                                                     
## (Intercept)                                         
## as.numeric(paste(redData$residual_rltl_tminus1)) ***
## redData$age_y                                    .  
## redData$feed_group2                                 
## redData$feed_groupNULL                              
## redData$genetic_groupS                              
## redData$birth_year                                  
## as.numeric(paste(redData$sample_interval))          
## redData$health_event_2week                          
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) a.((D$__ rdDt$g_ rdD$_2 rD$_NU rdD$_S rdDt$b_ a.((D$_)
## a.((D$__1)) -0.099                                                       
## redData$g_y -0.223  0.268                                                
## rdDt$fd_gr2 -0.015 -0.037   -0.068                                       
## rdDt$f_NULL  0.152 -0.055    0.006   0.392                               
## rdDt$gntc_S  0.115 -0.074   -0.012  -0.025  0.239                        
## rdDt$brth_y -1.000  0.099    0.223   0.015 -0.153 -0.115                 
## as.n((D$_)) -0.006 -0.010   -0.184   0.040 -0.027  0.053  0.004          
## rdDt$hlt__2 -0.003  0.045   -0.129   0.091  0.087  0.002  0.003   0.093
anova(mod1)
## Type III Analysis of Variance Table with Satterthwaite's method
##                                                   Sum Sq Mean Sq NumDF  DenDF
## as.numeric(paste(redData$residual_rltl_tminus1)) 0.44706 0.44706     1 978.44
## redData$age_y                                    0.04802 0.04802     1 680.44
## redData$feed_group                               0.00606 0.00303     2 161.75
## redData$genetic_group                            0.00852 0.00852     1 153.06
## redData$birth_year                               0.03073 0.03073     1 207.75
## as.numeric(paste(redData$sample_interval))       0.00173 0.00173     1 878.98
## redData$health_event_2week                       0.00090 0.00090     1 857.55
##                                                  F value    Pr(>F)    
## as.numeric(paste(redData$residual_rltl_tminus1)) 29.1317 8.486e-08 ***
## redData$age_y                                     3.1294   0.07734 .  
## redData$feed_group                                0.1975   0.82099    
## redData$genetic_group                             0.5553   0.45732    
## redData$birth_year                                2.0022   0.15857    
## as.numeric(paste(redData$sample_interval))        0.1126   0.73733    
## redData$health_event_2week                        0.0587   0.80862    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
mod1<-lmer(as.numeric(paste(redData$residual_rltl_t))~ as.numeric(paste(redData$residual_rltl_tminus1))+redData$age_y +   (1|redData$recoded_id), na.action=na.exclude)

summary(mod1)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: 
## as.numeric(paste(redData$residual_rltl_t)) ~ as.numeric(paste(redData$residual_rltl_tminus1)) +  
##     redData$age_y + (1 | redData$recoded_id)
## 
## REML criterion at convergence: -1114.7
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.47562 -0.65780 -0.07042  0.55074  3.06859 
## 
## Random effects:
##  Groups             Name        Variance Std.Dev.
##  redData$recoded_id (Intercept) 0.005188 0.07203 
##  Residual                       0.015458 0.12433 
## Number of obs: 1020, groups:  redData$recoded_id, 305
## 
## Fixed effects:
##                                                    Estimate Std. Error
## (Intercept)                                       -0.047332   0.009727
## as.numeric(paste(redData$residual_rltl_tminus1))   0.172457   0.028866
## redData$age_y                                      0.007691   0.003290
##                                                          df t value Pr(>|t|)
## (Intercept)                                      887.019161  -4.866 1.35e-06
## as.numeric(paste(redData$residual_rltl_tminus1)) 958.118312   5.974 3.25e-09
## redData$age_y                                    797.726818   2.337   0.0197
##                                                     
## (Intercept)                                      ***
## as.numeric(paste(redData$residual_rltl_tminus1)) ***
## redData$age_y                                    *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) a.((D$
## a.((D$__1)) -0.233       
## redData$g_y -0.804  0.264
anova(mod1)
## Type III Analysis of Variance Table with Satterthwaite's method
##                                                   Sum Sq Mean Sq NumDF  DenDF
## as.numeric(paste(redData$residual_rltl_tminus1)) 0.55177 0.55177     1 958.12
## redData$age_y                                    0.08446 0.08446     1 797.73
##                                                  F value    Pr(>F)    
## as.numeric(paste(redData$residual_rltl_tminus1)) 35.6935 3.253e-09 ***
## redData$age_y                                     5.4638   0.01966 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
mod1out<-as.data.frame(VarCorr(mod1),comp=c("Variance", "Std.Dev."))
repET<-mod1out[1,4]/sum(mod1out[1:nrow(mod1out),4])
repET
## [1] 0.2512912
#ggplot(redData, aes(x=age_d, y=as.numeric(paste(redData$diff_rltl_res))))+geom_point()+theme_gray(20)+ xlab("Age in days")+ ylab("RLTL change")

agePlot<- ggplot(redData, aes(x=as.factor(redData$age_y), y=as.numeric(paste(redData$diff_rltl_res))))+geom_jitter()+theme_gray(20)+ xlab("Age in years")+ ylab("RLTL change")+ geom_boxplot()

grid.arrange(SubsPlot, agePlot, ncol=2)
## Warning: Use of `redData$age_y` is discouraged. Use `age_y` instead.
## Warning: Use of `redData$diff_rltl_res` is discouraged. Use `diff_rltl_res`
## instead.
## Warning: Use of `redData$age_y` is discouraged. Use `age_y` instead.
## Warning: Use of `redData$diff_rltl_res` is discouraged. Use `diff_rltl_res`
## instead.

Investigate the effect of average telomere length over life on productive life span:

produc <-as.numeric(paste(uDeadData$amp_305_lactation))
## Warning: NAs introduced by coercion
mod2<-glm(as.numeric(paste(uDeadData$herd_life))~uDeadData$average_resid + produc)
summary(mod2)
## 
## Call:
## glm(formula = as.numeric(paste(uDeadData$herd_life)) ~ uDeadData$average_resid + 
##     produc)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -794.08  -361.00   -35.95   311.94  1283.75  
## 
## Coefficients:
##                           Estimate Std. Error t value Pr(>|t|)    
## (Intercept)             1205.21924   76.16160  15.825  < 2e-16 ***
## uDeadData$average_resid -135.44248  299.38143  -0.452    0.651    
## produc                     0.07095    0.01075   6.597 4.18e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 215047.8)
## 
##     Null deviance: 49655637  on 189  degrees of freedom
## Residual deviance: 40213930  on 187  degrees of freedom
##   (51 observations deleted due to missingness)
## AIC: 2877.1
## 
## Number of Fisher Scoring iterations: 2

Average telomere length over life is not significantly correlated with productive lifespan.

Average telomere change over life is associated with productive lifespan.

Repeat analysis for animals with more than two samples: uDeadDataMore2

redMore2<-read.csv("/Users/luise/Documents/Work/CattleTelomereLength/LongTelomereChangePaper/git/TelomereChangeInDairyCattle/more_2_sampl.csv")
u_red_more_2<-redMore2[!duplicated(redMore2$recoded_id),]

redDeadMore2<-subset(redMore2, redMore2$herd_life != "NULL")
uDeadDataMore2<-redDeadMore2[!duplicated(redDeadMore2$recoded_id),]


product<-as.numeric(paste(uDeadDataMore2$amp_305_lactation))
## Warning: NAs introduced by coercion
mod2<-glm(as.numeric(paste(uDeadDataMore2$herd_life))~uDeadDataMore2$average_resid+product)
summary(mod2)
## 
## Call:
## glm(formula = as.numeric(paste(uDeadDataMore2$herd_life)) ~ uDeadDataMore2$average_resid + 
##     product)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -774.53  -377.92   -42.55   335.50  1238.72  
## 
## Coefficients:
##                               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                  1.271e+03  8.264e+01  15.382  < 2e-16 ***
## uDeadDataMore2$average_resid 3.831e+00  3.035e+02   0.013     0.99    
## product                      6.308e-02  1.133e-02   5.566 9.25e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 215215.5)
## 
##     Null deviance: 45657680  on 183  degrees of freedom
## Residual deviance: 38954006  on 181  degrees of freedom
##   (29 observations deleted due to missingness)
## AIC: 2786.6
## 
## Number of Fisher Scoring iterations: 2
mod3a<-glm(as.numeric(paste(uDeadDataMore2$herd_life))~ as.numeric(paste(uDeadDataMore2$rltl_change))+product)

summary(mod3a)
## 
## Call:
## glm(formula = as.numeric(paste(uDeadDataMore2$herd_life)) ~ as.numeric(paste(uDeadDataMore2$rltl_change)) + 
##     product)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -787.33  -364.01   -52.38   318.81  1224.01  
## 
## Coefficients:
##                                                Estimate Std. Error t value
## (Intercept)                                   1.294e+03  8.446e+01  15.316
## as.numeric(paste(uDeadDataMore2$rltl_change)) 5.775e+02  6.028e+02   0.958
## product                                       6.238e-02  1.130e-02   5.521
##                                               Pr(>|t|)    
## (Intercept)                                    < 2e-16 ***
## as.numeric(paste(uDeadDataMore2$rltl_change))    0.339    
## product                                       1.16e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 214130)
## 
##     Null deviance: 45657680  on 183  degrees of freedom
## Residual deviance: 38757530  on 181  degrees of freedom
##   (29 observations deleted due to missingness)
## AIC: 2785.6
## 
## Number of Fisher Scoring iterations: 2
AIC(mod3a)
## [1] 2785.623
mod3b<-glm(as.numeric(paste(uDeadDataMore2$herd_life))~ as.numeric(paste(uDeadDataMore2$abs_rltl_change))+ product)
summary(mod3b)           
## 
## Call:
## glm(formula = as.numeric(paste(uDeadDataMore2$herd_life)) ~ as.numeric(paste(uDeadDataMore2$abs_rltl_change)) + 
##     product)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -810.56  -386.38   -32.43   326.01  1251.06  
## 
## Coefficients:
##                                                     Estimate Std. Error t value
## (Intercept)                                       1344.28343  110.91823  12.120
## as.numeric(paste(uDeadDataMore2$abs_rltl_change)) -541.34365  558.66779  -0.969
## product                                              0.06374    0.01129   5.643
##                                                   Pr(>|t|)    
## (Intercept)                                        < 2e-16 ***
## as.numeric(paste(uDeadDataMore2$abs_rltl_change))    0.334    
## product                                           6.34e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 214105)
## 
##     Null deviance: 45657680  on 183  degrees of freedom
## Residual deviance: 38753008  on 181  degrees of freedom
##   (29 observations deleted due to missingness)
## AIC: 2785.6
## 
## Number of Fisher Scoring iterations: 2
AIC(mod3b)
## [1] 2785.601
AIC(mod3b)-AIC(mod3a)
## [1] -0.02147028
#All together
mod3c<-glm(as.numeric(paste(uDeadDataMore2$herd_life))~ as.numeric(paste(uDeadDataMore2$rltl_change))+ as.numeric(paste(uDeadDataMore2$abs_rltl_change))+ uDeadDataMore2$average_resid +product)
             

summary(mod3c)
## 
## Call:
## glm(formula = as.numeric(paste(uDeadDataMore2$herd_life)) ~ as.numeric(paste(uDeadDataMore2$rltl_change)) + 
##     as.numeric(paste(uDeadDataMore2$abs_rltl_change)) + uDeadDataMore2$average_resid + 
##     product)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -813.50  -364.81   -61.04   328.87  1236.41  
## 
## Coefficients:
##                                                     Estimate Std. Error t value
## (Intercept)                                       1345.44375  114.08909  11.793
## as.numeric(paste(uDeadDataMore2$rltl_change))      443.74881  644.06446   0.689
## as.numeric(paste(uDeadDataMore2$abs_rltl_change)) -421.00056  595.89616  -0.706
## uDeadDataMore2$average_resid                         1.53645  311.59944   0.005
## product                                              0.06304    0.01140   5.532
##                                                   Pr(>|t|)    
## (Intercept)                                        < 2e-16 ***
## as.numeric(paste(uDeadDataMore2$rltl_change))        0.492    
## as.numeric(paste(uDeadDataMore2$abs_rltl_change))    0.481    
## uDeadDataMore2$average_resid                         0.996    
## product                                           1.11e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 215903.2)
## 
##     Null deviance: 45657680  on 183  degrees of freedom
## Residual deviance: 38646678  on 179  degrees of freedom
##   (29 observations deleted due to missingness)
## AIC: 2789.1
## 
## Number of Fisher Scoring iterations: 2
mod3d<-glm(as.numeric(paste(uDeadDataMore2$herd_life))~ as.numeric(paste(uDeadDataMore2$rltl_change))+ as.numeric(paste(uDeadDataMore2$abs_rltl_change))+ product)
             

summary(mod3d)
## 
## Call:
## glm(formula = as.numeric(paste(uDeadDataMore2$herd_life)) ~ as.numeric(paste(uDeadDataMore2$rltl_change)) + 
##     as.numeric(paste(uDeadDataMore2$abs_rltl_change)) + product)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -813.44  -364.72   -60.94   328.88  1236.64  
## 
## Coefficients:
##                                                     Estimate Std. Error t value
## (Intercept)                                       1345.32218  111.08303  12.111
## as.numeric(paste(uDeadDataMore2$rltl_change))      444.33049  631.40750   0.704
## as.numeric(paste(uDeadDataMore2$abs_rltl_change)) -420.49044  585.21439  -0.719
## product                                              0.06305    0.01135   5.554
##                                                   Pr(>|t|)    
## (Intercept)                                        < 2e-16 ***
## as.numeric(paste(uDeadDataMore2$rltl_change))        0.483    
## as.numeric(paste(uDeadDataMore2$abs_rltl_change))    0.473    
## product                                           9.92e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 214703.8)
## 
##     Null deviance: 45657680  on 183  degrees of freedom
## Residual deviance: 38646683  on 180  degrees of freedom
##   (29 observations deleted due to missingness)
## AIC: 2787.1
## 
## Number of Fisher Scoring iterations: 2

Cox proportional hazard model of productive lifespan with average lifetime RLTL as predictor

coxFit<-coxph(Surv(coxData$TimeMeasure, coxData$Event == 1)~coxData$average_resid)

summary(coxFit)
## Call:
## coxph(formula = Surv(coxData$TimeMeasure, coxData$Event == 1) ~ 
##     coxData$average_resid)
## 
##   n= 305, number of events= 241 
## 
##                         coef exp(coef) se(coef)     z Pr(>|z|)
## coxData$average_resid 0.3407    1.4059   0.5905 0.577    0.564
## 
##                       exp(coef) exp(-coef) lower .95 upper .95
## coxData$average_resid     1.406     0.7113    0.4419     4.473
## 
## Concordance= 0.509  (se = 0.021 )
## Likelihood ratio test= 0.33  on 1 df,   p=0.6
## Wald test            = 0.33  on 1 df,   p=0.6
## Score (logrank) test = 0.33  on 1 df,   p=0.6
summary(coxFit)$logtest["pvalue"]
##    pvalue 
## 0.5641153
Change to discrete scale to enable plotting
### discrete scale

coxFit<-coxph(Surv(coxData$TimeMeasure, coxData$Event == 1)~coxData$average_resid_3G)
summary(coxFit)
## Call:
## coxph(formula = Surv(coxData$TimeMeasure, coxData$Event == 1) ~ 
##     coxData$average_resid_3G)
## 
##   n= 305, number of events= 241 
## 
##                             coef exp(coef) se(coef)     z Pr(>|z|)
## coxData$average_resid_3G 0.01410   1.01420  0.07886 0.179    0.858
## 
##                          exp(coef) exp(-coef) lower .95 upper .95
## coxData$average_resid_3G     1.014      0.986     0.869     1.184
## 
## Concordance= 0.501  (se = 0.02 )
## Likelihood ratio test= 0.03  on 1 df,   p=0.9
## Wald test            = 0.03  on 1 df,   p=0.9
## Score (logrank) test = 0.03  on 1 df,   p=0.9
#ggsave(file = "/Users/luise/Documents/Work/CattleTelomereLength/AgingCell/Images/MeanRLTL.pdf", print(g2))

#write.csv(coxData, "/Users/luise/Documents/Work/CattleTelomereLength/AgingCell/rawdata/CaxData.csv")

Simple correlations of herd life with lifetime RLTL measures

They show the same as Cox proportional hazard models.

dData<-subset(data, data$herd_life != "NULL")
cor.test(as.numeric(paste(dData$herd_life)), dData$abs_rltl_change)
## 
##  Pearson's product-moment correlation
## 
## data:  as.numeric(paste(dData$herd_life)) and dData$abs_rltl_change
## t = -5.5853, df = 1020, p-value = 2.991e-08
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.2311491 -0.1121316
## sample estimates:
##       cor 
## -0.172269
cor.test(as.numeric(paste(dData$herd_life)), dData$rltl_change)
## 
##  Pearson's product-moment correlation
## 
## data:  as.numeric(paste(dData$herd_life)) and dData$rltl_change
## t = 8.4519, df = 1020, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.1976120 0.3122569
## sample estimates:
##       cor 
## 0.2558338
cor.test(as.numeric(paste(dData$herd_life)), dData$average_resid)
## 
##  Pearson's product-moment correlation
## 
## data:  as.numeric(paste(dData$herd_life)) and dData$average_resid
## t = -0.19531, df = 1020, p-value = 0.8452
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.06741195  0.05522731
## sample estimates:
##          cor 
## -0.006115315
### discrete scale
coxFit<-coxph(Surv(coxData$TimeMeasure, coxData$Event == 1)~coxData$abs_rltl_change_3G)

summary(coxFit)
## Call:
## coxph(formula = Surv(coxData$TimeMeasure, coxData$Event == 1) ~ 
##     coxData$abs_rltl_change_3G)
## 
##   n= 305, number of events= 241 
## 
##                               coef exp(coef) se(coef)    z Pr(>|z|)  
## coxData$abs_rltl_change_3G 0.17936   1.19645  0.08189 2.19   0.0285 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                            exp(coef) exp(-coef) lower .95 upper .95
## coxData$abs_rltl_change_3G     1.196     0.8358     1.019     1.405
## 
## Concordance= 0.542  (se = 0.02 )
## Likelihood ratio test= 4.8  on 1 df,   p=0.03
## Wald test            = 4.8  on 1 df,   p=0.03
## Score (logrank) test = 4.82  on 1 df,   p=0.03
coxFit<-survfit(Surv(coxData$TimeMeasure, coxData$Event == 1)~coxData$abs_rltl_change_3G)

summary(coxFit)
## Call: survfit(formula = Surv(coxData$TimeMeasure, coxData$Event == 
##     1) ~ coxData$abs_rltl_change_3G)
## 
##                 coxData$abs_rltl_change_3G=1 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##   495    102       1    0.990 0.00976       0.9713        1.000
##   501    101       1    0.980 0.01373       0.9539        1.000
##   636    100       1    0.971 0.01673       0.9383        1.000
##   641     99       1    0.961 0.01922       0.9238        0.999
##   648     98       1    0.951 0.02138       0.9100        0.994
##   705     97       1    0.941 0.02330       0.8966        0.988
##   724     96       1    0.931 0.02503       0.8836        0.982
##   727     95       1    0.922 0.02662       0.8708        0.975
##   735     94       1    0.912 0.02808       0.8583        0.969
##   794     93       1    0.902 0.02944       0.8461        0.962
##   805     92       1    0.892 0.03071       0.8339        0.954
##   830     91       1    0.882 0.03190       0.8220        0.947
##   831     90       1    0.873 0.03302       0.8102        0.940
##   833     89       1    0.863 0.03407       0.7985        0.932
##   836     88       1    0.853 0.03507       0.7869        0.925
##   895     87       1    0.843 0.03601       0.7754        0.917
##   897     86       1    0.833 0.03690       0.7641        0.909
##   904     85       1    0.824 0.03775       0.7528        0.901
##   909     84       1    0.814 0.03855       0.7416        0.893
##   956     83       1    0.804 0.03931       0.7304        0.885
##   995     82       1    0.794 0.04004       0.7194        0.877
##  1033     81       1    0.784 0.04072       0.7084        0.868
##  1082     80       1    0.775 0.04138       0.6975        0.860
##  1121     79       1    0.765 0.04200       0.6867        0.852
##  1128     78       1    0.755 0.04259       0.6759        0.843
##  1137     77       1    0.745 0.04315       0.6651        0.835
##  1158     76       2    0.725 0.04419       0.6439        0.817
##  1186     74       1    0.716 0.04466       0.6333        0.809
##  1209     73       1    0.706 0.04512       0.6228        0.800
##  1341     72       1    0.696 0.04554       0.6123        0.791
##  1356     71       1    0.686 0.04594       0.6019        0.782
##  1381     70       1    0.676 0.04632       0.5915        0.774
##  1483     69       1    0.667 0.04668       0.5812        0.765
##  1515     68       1    0.657 0.04701       0.5709        0.756
##  1551     67       1    0.647 0.04732       0.5607        0.747
##  1567     66       1    0.637 0.04761       0.5505        0.738
##  1573     65       1    0.627 0.04787       0.5403        0.729
##  1588     64       1    0.618 0.04812       0.5302        0.720
##  1596     63       1    0.608 0.04834       0.5201        0.710
##  1604     62       1    0.598 0.04855       0.5101        0.701
##  1612     61       1    0.588 0.04873       0.5001        0.692
##  1620     60       1    0.578 0.04889       0.4901        0.683
##  1624     59       1    0.569 0.04904       0.4802        0.673
##  1627     58       1    0.559 0.04916       0.4703        0.664
##  1633     57       1    0.549 0.04927       0.4605        0.655
##  1657     56       1    0.539 0.04935       0.4507        0.645
##  1766     55       1    0.529 0.04942       0.4409        0.636
##  1815     54       1    0.520 0.04947       0.4312        0.626
##  1825     53       1    0.510 0.04950       0.4215        0.617
##  1830     52       1    0.500 0.04951       0.4118        0.607
##  1834     51       1    0.490 0.04950       0.4022        0.597
##  1902     50       1    0.480 0.04947       0.3926        0.588
##  1903     49       1    0.471 0.04942       0.3830        0.578
##  1942     48       1    0.461 0.04935       0.3735        0.568
##  1960     47       1    0.451 0.04927       0.3641        0.559
##  2015     43       1    0.440 0.04923       0.3538        0.548
##  2017     42       1    0.430 0.04916       0.3437        0.538
##  2060     40       1    0.419 0.04909       0.3333        0.527
##  2072     39       1    0.409 0.04900       0.3229        0.517
##  2077     37       1    0.397 0.04890       0.3123        0.506
##  2140     35       1    0.386 0.04880       0.3014        0.495
##  2146     33       1    0.374 0.04871       0.2901        0.483
##  2177     31       1    0.362 0.04861       0.2786        0.471
##  2192     29       1    0.350 0.04851       0.2666        0.459
##  2214     28       1    0.337 0.04836       0.2547        0.447
##  2264     27       1    0.325 0.04816       0.2429        0.434
##  2275     26       1    0.312 0.04790       0.2313        0.422
##  2375     19       1    0.296 0.04812       0.2152        0.407
##  2437     15       1    0.276 0.04879       0.1954        0.390
##  2452     13       1    0.255 0.04944       0.1743        0.373
##  2457     12       1    0.234 0.04968       0.1541        0.354
##  2476     11       1    0.212 0.04950       0.1346        0.335
##  2567      7       1    0.182 0.05089       0.1053        0.315
##  2627      5       1    0.146 0.05214       0.0722        0.294
##  2708      4       1    0.109 0.05024       0.0444        0.269
## 
##                 coxData$abs_rltl_change_3G=2 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##   383    102       1   0.9902 0.00976       0.9713        1.000
##   499    101       1   0.9804 0.01373       0.9539        1.000
##   628    100       1   0.9706 0.01673       0.9383        1.000
##   658     99       1   0.9608 0.01922       0.9238        0.999
##   721     98       1   0.9510 0.02138       0.9100        0.994
##   731     97       1   0.9412 0.02330       0.8966        0.988
##   749     96       1   0.9314 0.02503       0.8836        0.982
##   822     95       1   0.9216 0.02662       0.8708        0.975
##   840     94       1   0.9118 0.02808       0.8583        0.969
##   908     93       1   0.9020 0.02944       0.8461        0.962
##   923     92       1   0.8922 0.03071       0.8339        0.954
##   928     91       1   0.8824 0.03190       0.8220        0.947
##   930     90       1   0.8725 0.03302       0.8102        0.940
##   948     89       1   0.8627 0.03407       0.7985        0.932
##   961     88       1   0.8529 0.03507       0.7869        0.925
##  1019     87       1   0.8431 0.03601       0.7754        0.917
##  1057     86       1   0.8333 0.03690       0.7641        0.909
##  1094     85       1   0.8235 0.03775       0.7528        0.901
##  1110     84       1   0.8137 0.03855       0.7416        0.893
##  1127     83       1   0.8039 0.03931       0.7304        0.885
##  1138     82       1   0.7941 0.04004       0.7194        0.877
##  1142     81       1   0.7843 0.04072       0.7084        0.868
##  1160     80       1   0.7745 0.04138       0.6975        0.860
##  1182     79       1   0.7647 0.04200       0.6867        0.852
##  1193     78       1   0.7549 0.04259       0.6759        0.843
##  1263     77       1   0.7451 0.04315       0.6651        0.835
##  1280     76       1   0.7353 0.04368       0.6545        0.826
##  1283     75       1   0.7255 0.04419       0.6439        0.817
##  1371     74       1   0.7157 0.04466       0.6333        0.809
##  1417     73       1   0.7059 0.04512       0.6228        0.800
##  1421     72       1   0.6961 0.04554       0.6123        0.791
##  1429     71       1   0.6863 0.04594       0.6019        0.782
##  1461     70       1   0.6765 0.04632       0.5915        0.774
##  1468     69       1   0.6667 0.04668       0.5812        0.765
##  1479     68       1   0.6569 0.04701       0.5709        0.756
##  1516     67       1   0.6471 0.04732       0.5607        0.747
##  1532     66       1   0.6373 0.04761       0.5505        0.738
##  1549     65       1   0.6275 0.04787       0.5403        0.729
##  1575     64       1   0.6176 0.04812       0.5302        0.720
##  1593     63       1   0.6078 0.04834       0.5201        0.710
##  1612     62       1   0.5980 0.04855       0.5101        0.701
##  1633     61       1   0.5882 0.04873       0.5001        0.692
##  1638     60       1   0.5784 0.04889       0.4901        0.683
##  1661     59       1   0.5686 0.04904       0.4802        0.673
##  1693     58       1   0.5588 0.04916       0.4703        0.664
##  1739     57       1   0.5490 0.04927       0.4605        0.655
##  1750     56       1   0.5392 0.04935       0.4507        0.645
##  1758     55       1   0.5294 0.04942       0.4409        0.636
##  1776     54       1   0.5196 0.04947       0.4312        0.626
##  1785     53       1   0.5098 0.04950       0.4215        0.617
##  1786     52       1   0.5000 0.04951       0.4118        0.607
##  1790     51       1   0.4902 0.04950       0.4022        0.597
##  1793     50       1   0.4804 0.04947       0.3926        0.588
##  1822     49       1   0.4706 0.04942       0.3830        0.578
##  1855     48       1   0.4608 0.04935       0.3735        0.568
##  1904     47       1   0.4510 0.04927       0.3641        0.559
##  1914     46       1   0.4412 0.04916       0.3546        0.549
##  1927     45       1   0.4314 0.04904       0.3452        0.539
##  1975     43       1   0.4213 0.04891       0.3356        0.529
##  1977     42       1   0.4113 0.04877       0.3260        0.519
##  1978     41       1   0.4013 0.04860       0.3165        0.509
##  1992     40       1   0.3912 0.04841       0.3070        0.499
##  2018     39       1   0.3812 0.04819       0.2975        0.488
##  2041     38       1   0.3712 0.04796       0.2881        0.478
##  2045     37       1   0.3611 0.04770       0.2788        0.468
##  2063     36       1   0.3511 0.04742       0.2695        0.458
##  2079     35       1   0.3411 0.04711       0.2602        0.447
##  2090     34       1   0.3311 0.04678       0.2510        0.437
##  2155     28       1   0.3192 0.04658       0.2398        0.425
##  2171     27       1   0.3074 0.04633       0.2288        0.413
##  2185     26       1   0.2956 0.04604       0.2178        0.401
##  2191     25       1   0.2838 0.04569       0.2070        0.389
##  2275     22       1   0.2709 0.04539       0.1950        0.376
##  2325     20       1   0.2573 0.04510       0.1825        0.363
##  2510     17       1   0.2422 0.04492       0.1684        0.348
##  2541     16       1   0.2270 0.04459       0.1545        0.334
##  2583     13       1   0.2096 0.04445       0.1383        0.318
##  2623     11       1   0.1905 0.04430       0.1208        0.301
##  2766      8       1   0.1667 0.04471       0.0986        0.282
##  2767      7       1   0.1429 0.04421       0.0779        0.262
##  2822      5       1   0.1143 0.04364       0.0541        0.242
##  2823      4       1   0.0857 0.04103       0.0336        0.219
## 
##                 coxData$abs_rltl_change_3G=3 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##    30    101       1   0.9901 0.00985       0.9710        1.000
##   273    100       1   0.9802 0.01386       0.9534        1.000
##   302     99       1   0.9703 0.01689       0.9377        1.000
##   330     98       1   0.9604 0.01941       0.9231        0.999
##   349     97       1   0.9505 0.02158       0.9091        0.994
##   476     96       1   0.9406 0.02352       0.8956        0.988
##   535     95       1   0.9307 0.02527       0.8825        0.982
##   544     94       1   0.9208 0.02687       0.8696        0.975
##   658     93       1   0.9109 0.02835       0.8570        0.968
##   699     92       1   0.9010 0.02972       0.8446        0.961
##   742     91       1   0.8911 0.03100       0.8324        0.954
##   750     90       1   0.8812 0.03220       0.8203        0.947
##   759     89       1   0.8713 0.03332       0.8084        0.939
##   760     88       1   0.8614 0.03438       0.7966        0.931
##   775     87       1   0.8515 0.03538       0.7849        0.924
##   798     86       1   0.8416 0.03633       0.7733        0.916
##   803     85       1   0.8317 0.03723       0.7618        0.908
##   845     84       1   0.8218 0.03808       0.7504        0.900
##   867     83       1   0.8119 0.03889       0.7391        0.892
##   888     82       1   0.8020 0.03965       0.7279        0.884
##   906     81       1   0.7921 0.04038       0.7168        0.875
##   919     80       1   0.7822 0.04107       0.7057        0.867
##   925     79       1   0.7723 0.04173       0.6947        0.859
##   932     78       1   0.7624 0.04235       0.6837        0.850
##   944     77       1   0.7525 0.04294       0.6728        0.842
##   958     76       1   0.7426 0.04350       0.6620        0.833
##   961     75       1   0.7327 0.04404       0.6513        0.824
##   982     74       1   0.7228 0.04454       0.6405        0.816
##  1006     73       1   0.7129 0.04502       0.6299        0.807
##  1045     72       1   0.7030 0.04547       0.6193        0.798
##  1069     71       1   0.6931 0.04589       0.6087        0.789
##  1142     70       1   0.6832 0.04629       0.5982        0.780
##  1147     69       1   0.6733 0.04667       0.5877        0.771
##  1156     68       1   0.6634 0.04702       0.5773        0.762
##  1165     67       1   0.6535 0.04735       0.5669        0.753
##  1177     66       1   0.6436 0.04766       0.5566        0.744
##  1195     65       1   0.6337 0.04794       0.5463        0.735
##  1205     64       1   0.6238 0.04820       0.5361        0.726
##  1222     63       1   0.6139 0.04844       0.5259        0.717
##  1301     62       1   0.6040 0.04866       0.5157        0.707
##  1345     61       1   0.5941 0.04886       0.5056        0.698
##  1364     60       1   0.5842 0.04904       0.4955        0.689
##  1373     59       1   0.5743 0.04920       0.4855        0.679
##  1378     58       1   0.5644 0.04934       0.4755        0.670
##  1384     57       1   0.5545 0.04946       0.4655        0.660
##  1436     56       1   0.5446 0.04955       0.4556        0.651
##  1530     55       1   0.5347 0.04963       0.4457        0.641
##  1557     54       1   0.5248 0.04969       0.4359        0.632
##  1561     53       1   0.5149 0.04973       0.4261        0.622
##  1565     52       1   0.5050 0.04975       0.4163        0.613
##  1597     51       1   0.4950 0.04975       0.4065        0.603
##  1608     50       1   0.4851 0.04973       0.3968        0.593
##  1609     49       1   0.4752 0.04969       0.3872        0.583
##  1610     48       1   0.4653 0.04963       0.3776        0.574
##  1615     47       1   0.4554 0.04955       0.3680        0.564
##  1671     45       1   0.4453 0.04948       0.3582        0.554
##  1701     44       1   0.4352 0.04938       0.3484        0.544
##  1715     43       1   0.4251 0.04925       0.3387        0.533
##  1756     42       1   0.4150 0.04911       0.3291        0.523
##  1762     41       1   0.4048 0.04894       0.3194        0.513
##  1766     40       1   0.3947 0.04876       0.3098        0.503
##  1770     39       1   0.3846 0.04854       0.3003        0.493
##  1786     38       1   0.3745 0.04831       0.2908        0.482
##  1848     37       1   0.3644 0.04805       0.2814        0.472
##  1904     36       1   0.3542 0.04777       0.2720        0.461
##  1937     35       1   0.3441 0.04747       0.2626        0.451
##  1940     34       1   0.3340 0.04714       0.2533        0.440
##  1994     33       1   0.3239 0.04678       0.2440        0.430
##  1999     32       1   0.3138 0.04640       0.2348        0.419
##  2001     31       1   0.3036 0.04600       0.2256        0.409
##  2035     27       1   0.2924 0.04565       0.2153        0.397
##  2060     26       1   0.2811 0.04526       0.2051        0.385
##  2067     25       2   0.2586 0.04434       0.1848        0.362
##  2075     23       1   0.2474 0.04382       0.1748        0.350
##  2086     22       1   0.2362 0.04324       0.1649        0.338
##  2113     18       1   0.2230 0.04279       0.1531        0.325
##  2166     16       1   0.2091 0.04232       0.1406        0.311
##  2235     14       1   0.1942 0.04185       0.1273        0.296
##  2285     11       1   0.1765 0.04160       0.1112        0.280
##  2333      8       1   0.1544 0.04185       0.0908        0.263
##  2486      4       1   0.1158 0.04586       0.0533        0.252
##  2579      3       1   0.0772 0.04392       0.0253        0.235
#### plot results for discrete scale

g3<- ggsurvplot(coxFit, data=coxData, pval = 0.029, ggtheme = theme_minimal(), risk.table.y.text.col = T, risk.table.y.text = FALSE, legend.labs=c("Tertile 1 (least absolute RLTL change)", "Tertile 2", "Tertile 3 (most absolute RLTL change)")) +theme_survminer(40)+ guides(colour = guide_legend(nrow = 3))
  
ggsurvplot(coxFit, data=coxData, pval = 0.029, ggtheme = theme_minimal(), risk.table.y.text.col = T, risk.table.y.text = FALSE, legend.labs=c("Tertile 1 (least abs RLTL change)", "Tertile 2", "Tertile 3 (most abs RLTL change)")) +theme_survminer(40)

#ggsave(file = "/Users/luise/Documents/Work/CattleTelomereLength/AgingCell/Images/AbsoluteRLTLchange.pdf", print(g3))
#continous scale

coxFit<-coxph(Surv(coxData$TimeMeasure, coxData$Event == 1)~coxData$rltl_change)

summary(coxFit)
## Call:
## coxph(formula = Surv(coxData$TimeMeasure, coxData$Event == 1) ~ 
##     coxData$rltl_change)
## 
##   n= 305, number of events= 241 
## 
##                          coef exp(coef)  se(coef)      z Pr(>|z|)    
## coxData$rltl_change -5.209383  0.005465  0.844821 -6.166 6.99e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                     exp(coef) exp(-coef) lower .95 upper .95
## coxData$rltl_change  0.005465        183  0.001043   0.02862
## 
## Concordance= 0.575  (se = 0.023 )
## Likelihood ratio test= 28.66  on 1 df,   p=9e-08
## Wald test            = 38.02  on 1 df,   p=7e-10
## Score (logrank) test = 33.59  on 1 df,   p=7e-09
#3 groups:
coxFit<-coxph(Surv(coxData$TimeMeasure, coxData$Event == 1)~coxData$rltl_change_3G)

summary(coxFit)
## Call:
## coxph(formula = Surv(coxData$TimeMeasure, coxData$Event == 1) ~ 
##     coxData$rltl_change_3G)
## 
##   n= 305, number of events= 241 
## 
##                            coef exp(coef) se(coef)      z Pr(>|z|)   
## coxData$rltl_change_3G -0.25745   0.77302  0.08656 -2.974  0.00294 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                        exp(coef) exp(-coef) lower .95 upper .95
## coxData$rltl_change_3G     0.773      1.294    0.6524    0.9159
## 
## Concordance= 0.56  (se = 0.019 )
## Likelihood ratio test= 8.88  on 1 df,   p=0.003
## Wald test            = 8.85  on 1 df,   p=0.003
## Score (logrank) test = 8.92  on 1 df,   p=0.003
coxFit<-survfit(Surv(coxData$TimeMeasure, coxData$Event == 1)~coxData$rltl_change_3G)

summary(coxFit)
## Call: survfit(formula = Surv(coxData$TimeMeasure, coxData$Event == 
##     1) ~ coxData$rltl_change_3G)
## 
##                 coxData$rltl_change_3G=1 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##    30    102       1   0.9902 0.00976       0.9713        1.000
##   273    101       1   0.9804 0.01373       0.9539        1.000
##   330    100       1   0.9706 0.01673       0.9383        1.000
##   349     99       1   0.9608 0.01922       0.9238        0.999
##   383     98       1   0.9510 0.02138       0.9100        0.994
##   476     97       1   0.9412 0.02330       0.8966        0.988
##   499     96       1   0.9314 0.02503       0.8836        0.982
##   501     95       1   0.9216 0.02662       0.8708        0.975
##   535     94       1   0.9118 0.02808       0.8583        0.969
##   544     93       1   0.9020 0.02944       0.8461        0.962
##   628     92       1   0.8922 0.03071       0.8339        0.954
##   636     91       1   0.8824 0.03190       0.8220        0.947
##   658     90       1   0.8725 0.03302       0.8102        0.940
##   699     89       1   0.8627 0.03407       0.7985        0.932
##   705     88       1   0.8529 0.03507       0.7869        0.925
##   724     87       1   0.8431 0.03601       0.7754        0.917
##   742     86       1   0.8333 0.03690       0.7641        0.909
##   750     85       1   0.8235 0.03775       0.7528        0.901
##   759     84       1   0.8137 0.03855       0.7416        0.893
##   760     83       1   0.8039 0.03931       0.7304        0.885
##   775     82       1   0.7941 0.04004       0.7194        0.877
##   845     81       1   0.7843 0.04072       0.7084        0.868
##   888     80       1   0.7745 0.04138       0.6975        0.860
##   906     79       1   0.7647 0.04200       0.6867        0.852
##   908     78       1   0.7549 0.04259       0.6759        0.843
##   923     77       1   0.7451 0.04315       0.6651        0.835
##   925     76       1   0.7353 0.04368       0.6545        0.826
##   944     75       1   0.7255 0.04419       0.6439        0.817
##   948     74       1   0.7157 0.04466       0.6333        0.809
##   961     73       2   0.6961 0.04554       0.6123        0.791
##   995     71       1   0.6863 0.04594       0.6019        0.782
##  1006     70       1   0.6765 0.04632       0.5915        0.774
##  1019     69       1   0.6667 0.04668       0.5812        0.765
##  1069     68       1   0.6569 0.04701       0.5709        0.756
##  1082     67       1   0.6471 0.04732       0.5607        0.747
##  1094     66       1   0.6373 0.04761       0.5505        0.738
##  1110     65       1   0.6275 0.04787       0.5403        0.729
##  1121     64       1   0.6176 0.04812       0.5302        0.720
##  1156     63       1   0.6078 0.04834       0.5201        0.710
##  1160     62       1   0.5980 0.04855       0.5101        0.701
##  1165     61       1   0.5882 0.04873       0.5001        0.692
##  1177     60       1   0.5784 0.04889       0.4901        0.683
##  1182     59       1   0.5686 0.04904       0.4802        0.673
##  1195     58       1   0.5588 0.04916       0.4703        0.664
##  1205     57       1   0.5490 0.04927       0.4605        0.655
##  1280     56       1   0.5392 0.04935       0.4507        0.645
##  1345     55       1   0.5294 0.04942       0.4409        0.636
##  1364     54       1   0.5196 0.04947       0.4312        0.626
##  1371     53       1   0.5098 0.04950       0.4215        0.617
##  1384     52       1   0.5000 0.04951       0.4118        0.607
##  1417     51       1   0.4902 0.04950       0.4022        0.597
##  1421     50       1   0.4804 0.04947       0.3926        0.588
##  1436     49       1   0.4706 0.04942       0.3830        0.578
##  1530     48       1   0.4608 0.04935       0.3735        0.568
##  1549     47       1   0.4510 0.04927       0.3641        0.559
##  1565     46       1   0.4412 0.04916       0.3546        0.549
##  1567     45       1   0.4314 0.04904       0.3452        0.539
##  1573     44       1   0.4216 0.04889       0.3358        0.529
##  1588     43       1   0.4118 0.04873       0.3265        0.519
##  1597     42       1   0.4020 0.04855       0.3172        0.509
##  1608     41       1   0.3922 0.04834       0.3080        0.499
##  1609     40       1   0.3824 0.04812       0.2988        0.489
##  1615     39       1   0.3725 0.04787       0.2896        0.479
##  1657     38       1   0.3627 0.04761       0.2805        0.469
##  1693     37       1   0.3529 0.04732       0.2714        0.459
##  1739     36       1   0.3431 0.04701       0.2623        0.449
##  1750     35       1   0.3333 0.04668       0.2533        0.439
##  1756     34       1   0.3235 0.04632       0.2444        0.428
##  1758     33       1   0.3137 0.04594       0.2354        0.418
##  1762     32       1   0.3039 0.04554       0.2266        0.408
##  1766     31       1   0.2941 0.04512       0.2177        0.397
##  1822     30       1   0.2843 0.04466       0.2090        0.387
##  1848     29       1   0.2745 0.04419       0.2002        0.376
##  1903     28       1   0.2647 0.04368       0.1916        0.366
##  1904     27       1   0.2549 0.04315       0.1829        0.355
##  1927     26       1   0.2451 0.04259       0.1744        0.345
##  1940     25       1   0.2353 0.04200       0.1658        0.334
##  1999     23       1   0.2251 0.04140       0.1569        0.323
##  2001     22       1   0.2148 0.04076       0.1481        0.312
##  2041     20       1   0.2041 0.04012       0.1388        0.300
##  2067     19       2   0.1826 0.03866       0.1206        0.277
##  2072     17       1   0.1719 0.03785       0.1116        0.265
##  2086     15       1   0.1604 0.03702       0.1020        0.252
##  2113     12       1   0.1470 0.03627       0.0907        0.238
##  2235      8       1   0.1287 0.03609       0.0742        0.223
##  2623      2       1   0.0643 0.04894       0.0145        0.286
##  2627      1       1   0.0000     NaN           NA           NA
## 
##                 coxData$rltl_change_3G=2 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##   641    102       1    0.990 0.00976       0.9713        1.000
##   727    101       1    0.980 0.01373       0.9539        1.000
##   731    100       1    0.971 0.01673       0.9383        1.000
##   749     99       1    0.961 0.01922       0.9238        0.999
##   794     98       1    0.951 0.02138       0.9100        0.994
##   805     97       1    0.941 0.02330       0.8966        0.988
##   904     96       1    0.931 0.02503       0.8836        0.982
##   909     95       1    0.922 0.02662       0.8708        0.975
##   928     94       1    0.912 0.02808       0.8583        0.969
##   956     93       1    0.902 0.02944       0.8461        0.962
##   958     92       1    0.892 0.03071       0.8339        0.954
##  1127     91       1    0.882 0.03190       0.8220        0.947
##  1137     90       1    0.873 0.03302       0.8102        0.940
##  1138     89       1    0.863 0.03407       0.7985        0.932
##  1142     88       2    0.843 0.03601       0.7754        0.917
##  1158     86       1    0.833 0.03690       0.7641        0.909
##  1222     85       1    0.824 0.03775       0.7528        0.901
##  1283     84       1    0.814 0.03855       0.7416        0.893
##  1341     83       1    0.804 0.03931       0.7304        0.885
##  1356     82       1    0.794 0.04004       0.7194        0.877
##  1378     81       1    0.784 0.04072       0.7084        0.868
##  1461     80       1    0.775 0.04138       0.6975        0.860
##  1468     79       1    0.765 0.04200       0.6867        0.852
##  1479     78       1    0.755 0.04259       0.6759        0.843
##  1532     77       1    0.745 0.04315       0.6651        0.835
##  1557     76       1    0.735 0.04368       0.6545        0.826
##  1596     75       1    0.725 0.04419       0.6439        0.817
##  1610     74       1    0.716 0.04466       0.6333        0.809
##  1612     73       1    0.706 0.04512       0.6228        0.800
##  1624     72       1    0.696 0.04554       0.6123        0.791
##  1627     71       1    0.686 0.04594       0.6019        0.782
##  1633     70       1    0.676 0.04632       0.5915        0.774
##  1661     69       1    0.667 0.04668       0.5812        0.765
##  1715     68       1    0.657 0.04701       0.5709        0.756
##  1766     67       1    0.647 0.04732       0.5607        0.747
##  1776     66       1    0.637 0.04761       0.5505        0.738
##  1786     65       1    0.627 0.04787       0.5403        0.729
##  1790     64       1    0.618 0.04812       0.5302        0.720
##  1793     63       1    0.608 0.04834       0.5201        0.710
##  1815     62       1    0.598 0.04855       0.5101        0.701
##  1825     61       1    0.588 0.04873       0.5001        0.692
##  1855     60       1    0.578 0.04889       0.4901        0.683
##  1937     59       1    0.569 0.04904       0.4802        0.673
##  1942     58       1    0.559 0.04916       0.4703        0.664
##  1977     57       1    0.549 0.04927       0.4605        0.655
##  1978     56       1    0.539 0.04935       0.4507        0.645
##  1992     55       1    0.529 0.04942       0.4409        0.636
##  2015     53       1    0.519 0.04949       0.4309        0.626
##  2018     52       1    0.509 0.04953       0.4210        0.616
##  2045     51       1    0.499 0.04956       0.4112        0.607
##  2060     50       1    0.489 0.04957       0.4013        0.597
##  2075     49       1    0.479 0.04955       0.3916        0.587
##  2079     48       1    0.469 0.04951       0.3818        0.577
##  2146     44       1    0.459 0.04953       0.3713        0.567
##  2166     43       1    0.448 0.04951       0.3609        0.556
##  2171     42       1    0.437 0.04947       0.3505        0.546
##  2177     41       1    0.427 0.04940       0.3402        0.535
##  2185     40       1    0.416 0.04930       0.3299        0.525
##  2191     38       1    0.405 0.04921       0.3194        0.514
##  2264     36       1    0.394 0.04911       0.3085        0.503
##  2275     33       1    0.382 0.04905       0.2970        0.491
##  2285     32       1    0.370 0.04895       0.2855        0.480
##  2375     26       1    0.356 0.04909       0.2715        0.466
##  2437     21       1    0.339 0.04959       0.2544        0.451
##  2457     19       1    0.321 0.05009       0.2365        0.436
##  2476     17       1    0.302 0.05057       0.2176        0.419
##  2486     16       1    0.283 0.05082       0.1993        0.403
##  2541     14       1    0.263 0.05106       0.1798        0.385
##  2567     11       1    0.239 0.05171       0.1565        0.365
##  2583     10       1    0.215 0.05178       0.1343        0.345
##  2766      7       1    0.184 0.05272       0.1053        0.323
##  2823      6       1    0.154 0.05213       0.0791        0.299
## 
##                 coxData$rltl_change_3G=3 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##   302    101       1   0.9901 0.00985       0.9710        1.000
##   495    100       1   0.9802 0.01386       0.9534        1.000
##   648     99       1   0.9703 0.01689       0.9377        1.000
##   658     98       1   0.9604 0.01941       0.9231        0.999
##   721     97       1   0.9505 0.02158       0.9091        0.994
##   735     96       1   0.9406 0.02352       0.8956        0.988
##   798     95       1   0.9307 0.02527       0.8825        0.982
##   803     94       1   0.9208 0.02687       0.8696        0.975
##   822     93       1   0.9109 0.02835       0.8570        0.968
##   830     92       1   0.9010 0.02972       0.8446        0.961
##   831     91       1   0.8911 0.03100       0.8324        0.954
##   833     90       1   0.8812 0.03220       0.8203        0.947
##   836     89       1   0.8713 0.03332       0.8084        0.939
##   840     88       1   0.8614 0.03438       0.7966        0.931
##   867     87       1   0.8515 0.03538       0.7849        0.924
##   895     86       1   0.8416 0.03633       0.7733        0.916
##   897     85       1   0.8317 0.03723       0.7618        0.908
##   919     84       1   0.8218 0.03808       0.7504        0.900
##   930     83       1   0.8119 0.03889       0.7391        0.892
##   932     82       1   0.8020 0.03965       0.7279        0.884
##   982     81       1   0.7921 0.04038       0.7168        0.875
##  1033     80       1   0.7822 0.04107       0.7057        0.867
##  1045     79       1   0.7723 0.04173       0.6947        0.859
##  1057     78       1   0.7624 0.04235       0.6837        0.850
##  1128     77       1   0.7525 0.04294       0.6728        0.842
##  1147     76       1   0.7426 0.04350       0.6620        0.833
##  1158     75       1   0.7327 0.04404       0.6513        0.824
##  1186     74       1   0.7228 0.04454       0.6405        0.816
##  1193     73       1   0.7129 0.04502       0.6299        0.807
##  1209     72       1   0.7030 0.04547       0.6193        0.798
##  1263     71       1   0.6931 0.04589       0.6087        0.789
##  1301     70       1   0.6832 0.04629       0.5982        0.780
##  1373     69       1   0.6733 0.04667       0.5877        0.771
##  1381     68       1   0.6634 0.04702       0.5773        0.762
##  1429     67       1   0.6535 0.04735       0.5669        0.753
##  1483     66       1   0.6436 0.04766       0.5566        0.744
##  1515     65       1   0.6337 0.04794       0.5463        0.735
##  1516     64       1   0.6238 0.04820       0.5361        0.726
##  1551     63       1   0.6139 0.04844       0.5259        0.717
##  1561     62       1   0.6040 0.04866       0.5157        0.707
##  1575     61       1   0.5941 0.04886       0.5056        0.698
##  1593     60       1   0.5842 0.04904       0.4955        0.689
##  1604     59       1   0.5743 0.04920       0.4855        0.679
##  1612     58       1   0.5644 0.04934       0.4755        0.670
##  1620     57       1   0.5545 0.04946       0.4655        0.660
##  1633     56       1   0.5446 0.04955       0.4556        0.651
##  1638     55       1   0.5347 0.04963       0.4457        0.641
##  1671     53       1   0.5246 0.04971       0.4356        0.632
##  1701     52       1   0.5145 0.04977       0.4256        0.622
##  1770     51       1   0.5044 0.04980       0.4156        0.612
##  1785     50       1   0.4943 0.04982       0.4057        0.602
##  1786     49       1   0.4842 0.04981       0.3958        0.592
##  1830     48       1   0.4741 0.04979       0.3859        0.582
##  1834     47       1   0.4640 0.04974       0.3761        0.573
##  1902     46       1   0.4540 0.04967       0.3663        0.563
##  1904     45       1   0.4439 0.04958       0.3566        0.552
##  1914     44       1   0.4338 0.04947       0.3469        0.542
##  1960     43       1   0.4237 0.04934       0.3372        0.532
##  1975     41       1   0.4134 0.04920       0.3273        0.522
##  1994     39       1   0.4028 0.04907       0.3172        0.511
##  2017     37       1   0.3919 0.04894       0.3068        0.501
##  2035     35       1   0.3807 0.04880       0.2961        0.489
##  2060     33       1   0.3691 0.04867       0.2851        0.478
##  2063     32       1   0.3576 0.04849       0.2741        0.466
##  2077     31       1   0.3461 0.04828       0.2633        0.455
##  2090     29       1   0.3341 0.04807       0.2520        0.443
##  2140     25       1   0.3208 0.04797       0.2393        0.430
##  2155     23       1   0.3068 0.04787       0.2260        0.417
##  2192     22       1   0.2929 0.04768       0.2129        0.403
##  2214     21       1   0.2789 0.04741       0.1999        0.389
##  2275     20       1   0.2650 0.04704       0.1871        0.375
##  2325     17       1   0.2494 0.04679       0.1727        0.360
##  2333     15       1   0.2328 0.04653       0.1573        0.344
##  2452     12       1   0.2134 0.04652       0.1392        0.327
##  2510     11       1   0.1940 0.04616       0.1217        0.309
##  2579      9       1   0.1724 0.04578       0.1025        0.290
##  2708      7       1   0.1478 0.04539       0.0810        0.270
##  2767      4       1   0.1108 0.04672       0.0485        0.253
##  2822      2       1   0.0554 0.04562       0.0110        0.278
#### plot results for discrete scale

g4<-ggsurvplot(coxFit, data=coxData, pval = 0.003, ggtheme = theme_minimal(), risk.table.y.text.col = T, risk.table.y.text = FALSE, legend.labs=c("Tertile 1 (most RLTL attrition)", "Tertile 2", "Tertile 3 (least RLTL attrition)"))+theme_survminer(40)+ guides(colour = guide_legend(nrow = 3))
  
  
ggsurvplot(coxFit, data=coxData, pval = 0.003, ggtheme = theme_minimal(), risk.table.y.text.col = T, risk.table.y.text = FALSE, legend.labs=c("Tertile 1 (most RLTL attrition)", "Tertile 2", "Tertile 3 (least RLTL attrition)")) +theme_survminer(40)

#ggsave(file = "/Users/luise/Documents/Work/CattleTelomereLength/AgingCell/Images/RLTLchange.pdf", print(g4))
ggplot(coxData, aes(x=as.factor(coxData$abs_rltl_change_3G), y= abs_rltl_change))+geom_boxplot() + xlab("Mean absolute RLTL change tertile") + ylab("Mean absolute RLTL change")+ theme_light(20)+ geom_jitter()
## Warning: Use of `coxData$abs_rltl_change_3G` is discouraged. Use
## `abs_rltl_change_3G` instead.

## Warning: Use of `coxData$abs_rltl_change_3G` is discouraged. Use
## `abs_rltl_change_3G` instead.

ggplot(coxData, aes(x=as.factor(coxData$rltl_change_3G), y= rltl_change))+geom_boxplot()+ xlab("Mean RLTL change tertile") + ylab("Mean RLTL change") + theme_light(20)+ geom_jitter()
## Warning: Use of `coxData$rltl_change_3G` is discouraged. Use `rltl_change_3G`
## instead.
## Warning: Use of `coxData$rltl_change_3G` is discouraged. Use `rltl_change_3G`
## instead.

ggplot(coxData, aes(x=as.factor(coxData$average_resid_3G), y= average_resid))+geom_boxplot()+ xlab("Mean RLTL tertile") + ylab("Mean RLTL")+ theme_light(20)+ geom_jitter()
## Warning: Use of `coxData$average_resid_3G` is discouraged. Use
## `average_resid_3G` instead.
## Warning: Use of `coxData$average_resid_3G` is discouraged. Use
## `average_resid_3G` instead.

RLTL at the age of 1 year and productive lifespan

coxFit<-coxph(Surv(coxData$TimeMeasure, coxData$Event == 1)~coxData$av_rltl_1_res)

summary(coxFit)
## Call:
## coxph(formula = Surv(coxData$TimeMeasure, coxData$Event == 1) ~ 
##     coxData$av_rltl_1_res)
## 
##   n= 284, number of events= 221 
##    (21 observations deleted due to missingness)
## 
##                          coef exp(coef) se(coef)      z Pr(>|z|)  
## coxData$av_rltl_1_res -0.9961    0.3693   0.4774 -2.087   0.0369 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                       exp(coef) exp(-coef) lower .95 upper .95
## coxData$av_rltl_1_res    0.3693      2.708    0.1449    0.9414
## 
## Concordance= 0.544  (se = 0.021 )
## Likelihood ratio test= 4.38  on 1 df,   p=0.04
## Wald test            = 4.35  on 1 df,   p=0.04
## Score (logrank) test = 4.36  on 1 df,   p=0.04

All three measures of lifetime RLTL dynamics and productive lifespan

coxFit<-coxph(Surv(coxData$TimeMeasure, coxData$Event == 1)~coxData$rltl_change+ coxData$average_resid+ coxData$abs_rltl_change)

summary(coxFit)

Including average milk production

coxFit<-coxph(Surv(coxDataR$TimeMeasure, coxDataR$Event == 1)~coxDataR$average_resid+ as.numeric(paste(coxDataR$amp_305_lactation)) + coxDataR$totalTELnumber)

summary(coxFit)
## Call:
## coxph(formula = Surv(coxDataR$TimeMeasure, coxDataR$Event == 
##     1) ~ coxDataR$average_resid + as.numeric(paste(coxDataR$amp_305_lactation)) + 
##     coxDataR$totalTELnumber)
## 
##   n= 253, number of events= 190 
## 
##                                                     coef  exp(coef)   se(coef)
## coxDataR$average_resid                         6.860e-01  1.986e+00  6.655e-01
## as.numeric(paste(coxDataR$amp_305_lactation)) -6.233e-05  9.999e-01  2.878e-05
## coxDataR$totalTELnumber                       -6.125e-01  5.420e-01  9.017e-02
##                                                    z Pr(>|z|)    
## coxDataR$average_resid                         1.031   0.3026    
## as.numeric(paste(coxDataR$amp_305_lactation)) -2.166   0.0303 *  
## coxDataR$totalTELnumber                       -6.793  1.1e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                                               exp(coef) exp(-coef) lower .95
## coxDataR$average_resid                           1.9857     0.5036    0.5388
## as.numeric(paste(coxDataR$amp_305_lactation))    0.9999     1.0001    0.9999
## coxDataR$totalTELnumber                          0.5420     1.8451    0.4542
##                                               upper .95
## coxDataR$average_resid                           7.3177
## as.numeric(paste(coxDataR$amp_305_lactation))    1.0000
## coxDataR$totalTELnumber                          0.6467
## 
## Concordance= 0.704  (se = 0.02 )
## Likelihood ratio test= 65.84  on 3 df,   p=3e-14
## Wald test            = 60.59  on 3 df,   p=4e-13
## Score (logrank) test = 62.99  on 3 df,   p=1e-13
coxFit<-coxph(Surv(coxDataR$TimeMeasure, coxDataR$Event == 1)~coxDataR$abs_rltl_change+ as.numeric(paste(coxDataR$amp_305_lactation))+coxDataR$totalTELnumber)

summary(coxFit)
## Call:
## coxph(formula = Surv(coxDataR$TimeMeasure, coxDataR$Event == 
##     1) ~ coxDataR$abs_rltl_change + as.numeric(paste(coxDataR$amp_305_lactation)) + 
##     coxDataR$totalTELnumber)
## 
##   n= 253, number of events= 190 
## 
##                                                     coef  exp(coef)   se(coef)
## coxDataR$abs_rltl_change                       1.066e+00  2.905e+00  1.147e+00
## as.numeric(paste(coxDataR$amp_305_lactation)) -6.366e-05  9.999e-01  2.868e-05
## coxDataR$totalTELnumber                       -5.948e-01  5.517e-01  9.143e-02
##                                                    z Pr(>|z|)    
## coxDataR$abs_rltl_change                       0.930   0.3524    
## as.numeric(paste(coxDataR$amp_305_lactation)) -2.219   0.0265 *  
## coxDataR$totalTELnumber                       -6.506 7.73e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                                               exp(coef) exp(-coef) lower .95
## coxDataR$abs_rltl_change                         2.9046     0.3443    0.3070
## as.numeric(paste(coxDataR$amp_305_lactation))    0.9999     1.0001    0.9999
## coxDataR$totalTELnumber                          0.5517     1.8127    0.4612
##                                               upper .95
## coxDataR$abs_rltl_change                        27.4796
## as.numeric(paste(coxDataR$amp_305_lactation))    1.0000
## coxDataR$totalTELnumber                          0.6599
## 
## Concordance= 0.703  (se = 0.02 )
## Likelihood ratio test= 65.62  on 3 df,   p=4e-14
## Wald test            = 60.45  on 3 df,   p=5e-13
## Score (logrank) test = 63.14  on 3 df,   p=1e-13
coxFit<-coxph(Surv(coxDataR$TimeMeasure, coxDataR$Event == 1)~coxDataR$rltl_change+ as.numeric(paste(coxDataR$amp_305_lactation))+coxDataR$totalTELnumber)

summary(coxFit)
## Call:
## coxph(formula = Surv(coxDataR$TimeMeasure, coxDataR$Event == 
##     1) ~ coxDataR$rltl_change + as.numeric(paste(coxDataR$amp_305_lactation)) + 
##     coxDataR$totalTELnumber)
## 
##   n= 253, number of events= 190 
## 
##                                                     coef  exp(coef)   se(coef)
## coxDataR$rltl_change                          -3.301e+00  3.685e-02  1.214e+00
## as.numeric(paste(coxDataR$amp_305_lactation)) -5.700e-05  9.999e-01  2.859e-05
## coxDataR$totalTELnumber                       -5.761e-01  5.621e-01  9.059e-02
##                                                    z Pr(>|z|)    
## coxDataR$rltl_change                          -2.719  0.00655 ** 
## as.numeric(paste(coxDataR$amp_305_lactation)) -1.994  0.04617 *  
## coxDataR$totalTELnumber                       -6.359 2.03e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                                               exp(coef) exp(-coef) lower .95
## coxDataR$rltl_change                            0.03685     27.140  0.003411
## as.numeric(paste(coxDataR$amp_305_lactation))   0.99994      1.000  0.999887
## coxDataR$totalTELnumber                         0.56211      1.779  0.470665
##                                               upper .95
## coxDataR$rltl_change                             0.3980
## as.numeric(paste(coxDataR$amp_305_lactation))    1.0000
## coxDataR$totalTELnumber                          0.6713
## 
## Concordance= 0.706  (se = 0.02 )
## Likelihood ratio test= 71.88  on 3 df,   p=2e-15
## Wald test            = 70.78  on 3 df,   p=3e-15
## Score (logrank) test = 72.11  on 3 df,   p=2e-15
coxFit<-coxph(Surv(coxDataR$TimeMeasure, coxDataR$Event == 1)~coxDataR$av_rltl_1_res + as.numeric(paste(coxDataR$amp_305_lactation))+ coxDataR$totalTELnumber)

summary(coxFit)
## Call:
## coxph(formula = Surv(coxDataR$TimeMeasure, coxDataR$Event == 
##     1) ~ coxDataR$av_rltl_1_res + as.numeric(paste(coxDataR$amp_305_lactation)) + 
##     coxDataR$totalTELnumber)
## 
##   n= 239, number of events= 177 
##    (14 observations deleted due to missingness)
## 
##                                                     coef  exp(coef)   se(coef)
## coxDataR$av_rltl_1_res                        -4.681e-02  9.543e-01  5.445e-01
## as.numeric(paste(coxDataR$amp_305_lactation)) -5.465e-05  9.999e-01  2.984e-05
## coxDataR$totalTELnumber                       -5.417e-01  5.818e-01  9.348e-02
##                                                    z Pr(>|z|)    
## coxDataR$av_rltl_1_res                        -0.086   0.9315    
## as.numeric(paste(coxDataR$amp_305_lactation)) -1.831   0.0671 .  
## coxDataR$totalTELnumber                       -5.794 6.86e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                                               exp(coef) exp(-coef) lower .95
## coxDataR$av_rltl_1_res                           0.9543      1.048    0.3282
## as.numeric(paste(coxDataR$amp_305_lactation))    0.9999      1.000    0.9999
## coxDataR$totalTELnumber                          0.5818      1.719    0.4844
##                                               upper .95
## coxDataR$av_rltl_1_res                           2.7745
## as.numeric(paste(coxDataR$amp_305_lactation))    1.0000
## coxDataR$totalTELnumber                          0.6988
## 
## Concordance= 0.686  (se = 0.021 )
## Likelihood ratio test= 48.63  on 3 df,   p=2e-10
## Wald test            = 44.57  on 3 df,   p=1e-09
## Score (logrank) test = 46.14  on 3 df,   p=5e-10
#total TEL sample number was not significant and therefore will be eliminated:

I THINK I SHOULD USE THE MODELS BELOW

coxFit<-coxph(Surv(coxDataR$TimeMeasure, coxDataR$Event == 1)~coxDataR$average_resid+ as.numeric(paste(coxDataR$amp_305_lactation)))

summary(coxFit)
## Call:
## coxph(formula = Surv(coxDataR$TimeMeasure, coxDataR$Event == 
##     1) ~ coxDataR$average_resid + as.numeric(paste(coxDataR$amp_305_lactation)))
## 
##   n= 253, number of events= 190 
## 
##                                                     coef  exp(coef)   se(coef)
## coxDataR$average_resid                         2.315e-01  1.260e+00  6.708e-01
## as.numeric(paste(coxDataR$amp_305_lactation)) -9.933e-05  9.999e-01  2.822e-05
##                                                    z Pr(>|z|)    
## coxDataR$average_resid                         0.345 0.730011    
## as.numeric(paste(coxDataR$amp_305_lactation)) -3.520 0.000432 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                                               exp(coef) exp(-coef) lower .95
## coxDataR$average_resid                           1.2605     0.7934    0.3385
## as.numeric(paste(coxDataR$amp_305_lactation))    0.9999     1.0001    0.9998
##                                               upper .95
## coxDataR$average_resid                            4.693
## as.numeric(paste(coxDataR$amp_305_lactation))     1.000
## 
## Concordance= 0.598  (se = 0.023 )
## Likelihood ratio test= 11.93  on 2 df,   p=0.003
## Wald test            = 12.43  on 2 df,   p=0.002
## Score (logrank) test = 12.44  on 2 df,   p=0.002
coxFit<-coxph(Surv(coxDataR$TimeMeasure, coxDataR$Event == 1)~coxDataR$abs_rltl_change+ as.numeric(paste(coxDataR$amp_305_lactation)))

summary(coxFit)
## Call:
## coxph(formula = Surv(coxDataR$TimeMeasure, coxDataR$Event == 
##     1) ~ coxDataR$abs_rltl_change + as.numeric(paste(coxDataR$amp_305_lactation)))
## 
##   n= 253, number of events= 190 
## 
##                                                     coef  exp(coef)   se(coef)
## coxDataR$abs_rltl_change                       2.411e+00  1.115e+01  1.146e+00
## as.numeric(paste(coxDataR$amp_305_lactation)) -1.018e-04  9.999e-01  2.818e-05
##                                                    z Pr(>|z|)    
## coxDataR$abs_rltl_change                       2.104 0.035391 *  
## as.numeric(paste(coxDataR$amp_305_lactation)) -3.611 0.000305 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                                               exp(coef) exp(-coef) lower .95
## coxDataR$abs_rltl_change                        11.1466    0.08971    1.1793
## as.numeric(paste(coxDataR$amp_305_lactation))    0.9999    1.00010    0.9998
##                                               upper .95
## coxDataR$abs_rltl_change                          105.4
## as.numeric(paste(coxDataR$amp_305_lactation))       1.0
## 
## Concordance= 0.598  (se = 0.023 )
## Likelihood ratio test= 15.99  on 2 df,   p=3e-04
## Wald test            = 17.06  on 2 df,   p=2e-04
## Score (logrank) test = 17.08  on 2 df,   p=2e-04
coxFit<-coxph(Surv(coxDataR$TimeMeasure, coxDataR$Event == 1)~coxDataR$rltl_change+ as.numeric(paste(coxDataR$amp_305_lactation)))

summary(coxFit)
## Call:
## coxph(formula = Surv(coxDataR$TimeMeasure, coxDataR$Event == 
##     1) ~ coxDataR$rltl_change + as.numeric(paste(coxDataR$amp_305_lactation)))
## 
##   n= 253, number of events= 190 
## 
##                                                     coef  exp(coef)   se(coef)
## coxDataR$rltl_change                          -5.035e+00  6.506e-03  1.315e+00
## as.numeric(paste(coxDataR$amp_305_lactation)) -9.057e-05  9.999e-01  2.819e-05
##                                                    z Pr(>|z|)    
## coxDataR$rltl_change                          -3.829 0.000129 ***
## as.numeric(paste(coxDataR$amp_305_lactation)) -3.212 0.001317 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                                               exp(coef) exp(-coef) lower .95
## coxDataR$rltl_change                           0.006506      153.7 0.0004944
## as.numeric(paste(coxDataR$amp_305_lactation))  0.999909        1.0 0.9998542
##                                               upper .95
## coxDataR$rltl_change                            0.08562
## as.numeric(paste(coxDataR$amp_305_lactation))   0.99996
## 
## Concordance= 0.601  (se = 0.022 )
## Likelihood ratio test= 25.19  on 2 df,   p=3e-06
## Wald test            = 28.18  on 2 df,   p=8e-07
## Score (logrank) test = 26.7  on 2 df,   p=2e-06
test RLTL at the age of 1 year in the same model as mean RLTL change
coxFit<-coxph(Surv(coxDataR$TimeMeasure, coxDataR$Event == 1)~coxDataR$rltl_change+ coxDataR$av_rltl_1_res+ as.numeric(paste(coxDataR$amp_305_lactation)))

summary(coxFit)
## Call:
## coxph(formula = Surv(coxDataR$TimeMeasure, coxDataR$Event == 
##     1) ~ coxDataR$rltl_change + coxDataR$av_rltl_1_res + as.numeric(paste(coxDataR$amp_305_lactation)))
## 
##   n= 239, number of events= 177 
##    (14 observations deleted due to missingness)
## 
##                                                     coef  exp(coef)   se(coef)
## coxDataR$rltl_change                          -5.030e+00  6.537e-03  1.379e+00
## coxDataR$av_rltl_1_res                        -6.552e-01  5.193e-01  5.263e-01
## as.numeric(paste(coxDataR$amp_305_lactation)) -7.896e-05  9.999e-01  2.928e-05
##                                                    z Pr(>|z|)    
## coxDataR$rltl_change                          -3.648 0.000264 ***
## coxDataR$av_rltl_1_res                        -1.245 0.213134    
## as.numeric(paste(coxDataR$amp_305_lactation)) -2.696 0.007008 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                                               exp(coef) exp(-coef) lower .95
## coxDataR$rltl_change                           0.006537    152.972 0.0004382
## coxDataR$av_rltl_1_res                         0.519327      1.926 0.1851255
## as.numeric(paste(coxDataR$amp_305_lactation))  0.999921      1.000 0.9998637
##                                               upper .95
## coxDataR$rltl_change                            0.09752
## coxDataR$av_rltl_1_res                          1.45685
## as.numeric(paste(coxDataR$amp_305_lactation))   0.99998
## 
## Concordance= 0.595  (se = 0.023 )
## Likelihood ratio test= 22.47  on 3 df,   p=5e-05
## Wald test            = 25.18  on 3 df,   p=1e-05
## Score (logrank) test = 23.85  on 3 df,   p=3e-05
All three measures of lifetime RLTL dynamics in the same model including milk productivity measure:
coxFit<-coxph(Surv(coxDataR$TimeMeasure, coxDataR$Event == 1)~coxDataR$rltl_change+coxDataR$average_resid+coxDataR$abs_rltl_change+ as.numeric(paste(coxDataR$amp_305_lactation)))

summary(coxFit)
## Call:
## coxph(formula = Surv(coxDataR$TimeMeasure, coxDataR$Event == 
##     1) ~ coxDataR$rltl_change + coxDataR$average_resid + coxDataR$abs_rltl_change + 
##     as.numeric(paste(coxDataR$amp_305_lactation)))
## 
##   n= 253, number of events= 190 
## 
##                                                     coef  exp(coef)   se(coef)
## coxDataR$rltl_change                          -4.525e+00  1.083e-02  1.387e+00
## coxDataR$average_resid                        -6.968e-02  9.327e-01  6.840e-01
## coxDataR$abs_rltl_change                       1.237e+00  3.447e+00  1.224e+00
## as.numeric(paste(coxDataR$amp_305_lactation)) -9.162e-05  9.999e-01  2.824e-05
##                                                    z Pr(>|z|)   
## coxDataR$rltl_change                          -3.262  0.00111 **
## coxDataR$average_resid                        -0.102  0.91885   
## coxDataR$abs_rltl_change                       1.011  0.31198   
## as.numeric(paste(coxDataR$amp_305_lactation)) -3.244  0.00118 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                                               exp(coef) exp(-coef) lower .95
## coxDataR$rltl_change                            0.01083    92.3027 0.0007142
## coxDataR$average_resid                          0.93269     1.0722 0.2440869
## coxDataR$abs_rltl_change                        3.44666     0.2901 0.3130755
## as.numeric(paste(coxDataR$amp_305_lactation))   0.99991     1.0001 0.9998530
##                                               upper .95
## coxDataR$rltl_change                             0.1643
## coxDataR$average_resid                           3.5639
## coxDataR$abs_rltl_change                        37.9445
## as.numeric(paste(coxDataR$amp_305_lactation))    1.0000
## 
## Concordance= 0.6  (se = 0.022 )
## Likelihood ratio test= 26.21  on 4 df,   p=3e-05
## Wald test            = 30.41  on 4 df,   p=4e-06
## Score (logrank) test = 28.88  on 4 df,   p=8e-06
Repeat Cox proportional hazard model without without the fist measurement.
red_data_l_new<-read.csv("/Users/luise/Documents/Work/CattleTelomereLength/LongTelomereChangePaper/git/TelomereChangeInDairyCattle/excluded_less_3_month.csv")
uniq_red_data_l_new<-red_data_l_new[!duplicated(red_data_l_new$recoded_id),]

Histogram of samples per animal after removing early measurements

hist1b<-ggplot(uniq_red_data_l_new, aes(x = uniq_red_data_l_new$SamplesPerAnimal)) + 
    geom_histogram(colour= "black", fill = mycoloursP[15], binwidth = 1) + xlab("Number of samples per animal") +  
  geom_vline(xintercept= mean(as.numeric(paste(uniq_red_data_l_new$SamplesPerAnimal)), na.rm=TRUE), colour="red")+ 
  annotate("text", x = 0.7, y = 100, label = "(A)", fontface =2, size = 10) + theme_classic(20)+
  scale_x_continuous(breaks=c(1:8))
## Warning in mean(as.numeric(paste(uniq_red_data_l_new$SamplesPerAnimal)), : NAs
## introduced by coercion
hist1b
## Warning: Use of `uniq_red_data_l_new$SamplesPerAnimal` is discouraged. Use
## `SamplesPerAnimal` instead.
## Warning: Removed 27 rows containing non-finite values (stat_bin).

For dataframe that excludes the first sample

hist2b<-ggplot(red_data_l_new, aes(x = red_data_l_new$residual_rltl_t)) + 
    geom_histogram(colour= "black", fill = mycoloursP[16]) + xlab("Adjusted RLTL") + 
  geom_vline(xintercept= mean(as.numeric(paste(red_data_l_new$residual_rltl_t))), colour="red") + 
  annotate("text", x = -0.45, y = 120, label = "(B)", fontface =2, size = 10)+
  theme_classic(20)+
    stat_function(
        fun = function(x, mean, sd, n){
             n/27* dnorm(x = x, mean = mean, sd = sd)
        }, 
        args = with(red_data_l_new, c(mean = mean(red_data_l_new$residual_rltl_t), sd = sd(red_data_l_new$residual_rltl_t), n
= length(red_data_l_new$residual_rltl_t)))
    ) 

hist2b
## Warning: Use of `red_data_l_new$residual_rltl_t` is discouraged. Use
## `residual_rltl_t` instead.

## Warning: Use of `red_data_l_new$residual_rltl_t` is discouraged. Use
## `residual_rltl_t` instead.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

#red_data_l_new_ch<-subset(red_data_l_new, red_data_l_new$SamplesPerAnimal != 2)


diff<- as.numeric(paste(red_data_l_new$diff_rltl_res))

hist3b<-ggplot(red_data_l_new, aes(x = as.numeric(paste(red_data_l_new$diff_rltl_res)))) + 
    geom_histogram(colour= "black", fill = mycoloursP[17]) + xlab("Adjusted RLTL change") +  
  geom_vline(xintercept= mean(as.numeric(paste(red_data_l_new$diff_rltl_res))), colour="red") +
  annotate("text", x = -0.75, y = 95, label = "(C)", fontface =2, size= 10)+
  theme_classic(20)+
    stat_function(
        fun = function(x, mean, sd, n){
             n/23* dnorm(x = x, mean = mean, sd = sd)
        }, 
        args = with(red_data_l_new, c(mean = mean(as.numeric(paste(red_data_l_new$diff_rltl_res))), sd = sd(as.numeric(paste(red_data_l_new$diff_rltl_res))), n
= length(as.numeric(paste(red_data_l_new$diff_rltl_res))))))

hist3b
## Warning: Use of `red_data_l_new$diff_rltl_res` is discouraged. Use
## `diff_rltl_res` instead.

## Warning: Use of `red_data_l_new$diff_rltl_res` is discouraged. Use
## `diff_rltl_res` instead.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

The plot above looks so much like hist1 that I think I may have made a mistake. when I look for summary output of diff_rltl_res in both the complete dataset (data) as well as the subsetted dataset (red_data_l_new) summary statistics are really very similar, but not identical. By removing single samples, I should have removed samples without a first change measurement. I have to think about this for a little bit.

summary(as.numeric(paste(red_data_l_new$diff_rltl_res)))
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## -0.72730 -0.15316 -0.03555 -0.03257  0.08741  0.48019
summary(as.numeric(paste(data$diff_rltl_res)))
## Warning in summary(as.numeric(paste(data$diff_rltl_res))): NAs introduced by
## coercion
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max.     NA's 
## -0.72730 -0.15310 -0.03538 -0.03243  0.08796  0.48019      305
u_red_data_l_new<-red_data_l_new[!duplicated(red_data_l_new$recoded_id),]


coxDataL<-data.frame(recoded_id=u_red_data_l_new$recoded_id, dob=u_red_data_l_new$dob, herd_life=u_red_data_l_new$herd_life, cull_date=u_red_data_l_new$cull_date, average_resid=u_red_data_l_new$average_resid, rltl_change=u_red_data_l_new$rltl_change, abs_rltl_change=u_red_data_l_new$abs_rltl_change, amp_305_lactation= u_red_data_l_new$amp_305_lactation)


coxDataL$CutoffDate<-"01/01/2017"
coxDataL$Event<-ifelse(coxDataL$cull_date=="NULL", 0, 1)
coxDataL$TimeMeasure<-ifelse(coxDataL$Event==0, as.Date(coxDataL$CutoffDate, "%d/%m/%Y")- as.Date(coxDataL$dob, "%d/%m/%Y"), as.numeric(paste(coxDataL$herd_life)))
## Warning in ifelse(coxDataL$Event == 0, as.Date(coxDataL$CutoffDate, "%d/%m/%Y")
## - : NAs introduced by coercion
coxFit<-coxph(Surv(coxDataL$TimeMeasure, coxDataL$Event == 1)~coxDataL$average_resid)

summary(coxFit)
## Call:
## coxph(formula = Surv(coxDataL$TimeMeasure, coxDataL$Event == 
##     1) ~ coxDataL$average_resid)
## 
##   n= 304, number of events= 240 
## 
##                           coef exp(coef) se(coef)      z Pr(>|z|)
## coxDataL$average_resid -0.5647    0.5686   0.5667 -0.996    0.319
## 
##                        exp(coef) exp(-coef) lower .95 upper .95
## coxDataL$average_resid    0.5686      1.759    0.1873     1.726
## 
## Concordance= 0.524  (se = 0.021 )
## Likelihood ratio test= 1  on 1 df,   p=0.3
## Wald test            = 0.99  on 1 df,   p=0.3
## Score (logrank) test = 0.99  on 1 df,   p=0.3
coxFit<-coxph(Surv(coxDataL$TimeMeasure, coxDataL$Event == 1)~coxDataL$abs_rltl_change)

summary(coxFit)
## Call:
## coxph(formula = Surv(coxDataL$TimeMeasure, coxDataL$Event == 
##     1) ~ coxDataL$abs_rltl_change)
## 
##   n= 304, number of events= 240 
## 
##                             coef exp(coef) se(coef)     z Pr(>|z|)   
## coxDataL$abs_rltl_change  2.7934   16.3372   0.9814 2.846  0.00442 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                          exp(coef) exp(-coef) lower .95 upper .95
## coxDataL$abs_rltl_change     16.34    0.06121     2.387     111.8
## 
## Concordance= 0.537  (se = 0.023 )
## Likelihood ratio test= 7.49  on 1 df,   p=0.006
## Wald test            = 8.1  on 1 df,   p=0.004
## Score (logrank) test = 8.05  on 1 df,   p=0.005
coxFit<-coxph(Surv(coxDataL$TimeMeasure, coxDataL$Event == 1)~coxDataL$rltl_change)

summary(coxFit)
## Call:
## coxph(formula = Surv(coxDataL$TimeMeasure, coxDataL$Event == 
##     1) ~ coxDataL$rltl_change)
## 
##   n= 304, number of events= 240 
## 
##                           coef exp(coef)  se(coef)      z Pr(>|z|)    
## coxDataL$rltl_change -5.065332  0.006312  0.869793 -5.824 5.76e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                      exp(coef) exp(-coef) lower .95 upper .95
## coxDataL$rltl_change  0.006312      158.4  0.001148   0.03472
## 
## Concordance= 0.572  (se = 0.023 )
## Likelihood ratio test= 26.01  on 1 df,   p=3e-07
## Wald test            = 33.91  on 1 df,   p=6e-09
## Score (logrank) test = 30.05  on 1 df,   p=4e-08

#now without first measurement but including milk productivity

#now I will include average milk production, TEL sample number was tested, but removed from the models because it was not statistically associated with productive lifespan

coxDataRL<-subset(coxDataL, coxDataL$amp_305_lactation != "NULL")

#coxDataRL<-merge(coxDataRL, maxSN, by = "recoded_id")


coxFit<-coxph(Surv(coxDataRL$TimeMeasure, coxDataRL$Event == 1)~coxDataRL$average_resid+ as.numeric(paste(coxDataRL$amp_305_lactation)))

summary(coxFit)
## Call:
## coxph(formula = Surv(coxDataRL$TimeMeasure, coxDataRL$Event == 
##     1) ~ coxDataRL$average_resid + as.numeric(paste(coxDataRL$amp_305_lactation)))
## 
##   n= 253, number of events= 190 
## 
##                                                      coef  exp(coef)   se(coef)
## coxDataRL$average_resid                        -4.847e-01  6.159e-01  6.494e-01
## as.numeric(paste(coxDataRL$amp_305_lactation)) -9.853e-05  9.999e-01  2.812e-05
##                                                     z Pr(>|z|)    
## coxDataRL$average_resid                        -0.746 0.455445    
## as.numeric(paste(coxDataRL$amp_305_lactation)) -3.504 0.000458 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                                                exp(coef) exp(-coef) lower .95
## coxDataRL$average_resid                           0.6159      1.624    0.1725
## as.numeric(paste(coxDataRL$amp_305_lactation))    0.9999      1.000    0.9998
##                                                upper .95
## coxDataRL$average_resid                            2.199
## as.numeric(paste(coxDataRL$amp_305_lactation))     1.000
## 
## Concordance= 0.599  (se = 0.023 )
## Likelihood ratio test= 12.37  on 2 df,   p=0.002
## Wald test            = 12.93  on 2 df,   p=0.002
## Score (logrank) test = 12.92  on 2 df,   p=0.002
coxFit<-coxph(Surv(coxDataRL$TimeMeasure, coxDataRL$Event == 1)~coxDataRL$abs_rltl_change+ as.numeric(paste(coxDataRL$amp_305_lactation)))

summary(coxFit)
## Call:
## coxph(formula = Surv(coxDataRL$TimeMeasure, coxDataRL$Event == 
##     1) ~ coxDataRL$abs_rltl_change + as.numeric(paste(coxDataRL$amp_305_lactation)))
## 
##   n= 253, number of events= 190 
## 
##                                                      coef  exp(coef)   se(coef)
## coxDataRL$abs_rltl_change                       2.403e+00  1.105e+01  1.147e+00
## as.numeric(paste(coxDataRL$amp_305_lactation)) -1.018e-04  9.999e-01  2.818e-05
##                                                     z Pr(>|z|)    
## coxDataRL$abs_rltl_change                       2.095 0.036134 *  
## as.numeric(paste(coxDataRL$amp_305_lactation)) -3.610 0.000306 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                                                exp(coef) exp(-coef) lower .95
## coxDataRL$abs_rltl_change                        11.0524    0.09048    1.1680
## as.numeric(paste(coxDataRL$amp_305_lactation))    0.9999    1.00010    0.9998
##                                                upper .95
## coxDataRL$abs_rltl_change                          104.6
## as.numeric(paste(coxDataRL$amp_305_lactation))       1.0
## 
## Concordance= 0.598  (se = 0.023 )
## Likelihood ratio test= 15.95  on 2 df,   p=3e-04
## Wald test            = 17.02  on 2 df,   p=2e-04
## Score (logrank) test = 17.04  on 2 df,   p=2e-04
coxFit<-coxph(Surv(coxDataRL$TimeMeasure, coxDataRL$Event == 1)~coxDataRL$rltl_change + as.numeric(paste(coxDataRL$amp_305_lactation)))

summary(coxFit)
## Call:
## coxph(formula = Surv(coxDataRL$TimeMeasure, coxDataRL$Event == 
##     1) ~ coxDataRL$rltl_change + as.numeric(paste(coxDataRL$amp_305_lactation)))
## 
##   n= 253, number of events= 190 
## 
##                                                      coef  exp(coef)   se(coef)
## coxDataRL$rltl_change                          -5.056e+00  6.370e-03  1.315e+00
## as.numeric(paste(coxDataRL$amp_305_lactation)) -9.096e-05  9.999e-01  2.818e-05
##                                                     z Pr(>|z|)    
## coxDataRL$rltl_change                          -3.844 0.000121 ***
## as.numeric(paste(coxDataRL$amp_305_lactation)) -3.228 0.001247 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                                                exp(coef) exp(-coef) lower .95
## coxDataRL$rltl_change                            0.00637        157 0.0004836
## as.numeric(paste(coxDataRL$amp_305_lactation))   0.99991          1 0.9998538
##                                                upper .95
## coxDataRL$rltl_change                             0.0839
## as.numeric(paste(coxDataRL$amp_305_lactation))    1.0000
## 
## Concordance= 0.601  (se = 0.022 )
## Likelihood ratio test= 25.3  on 2 df,   p=3e-06
## Wald test            = 28.32  on 2 df,   p=7e-07
## Score (logrank) test = 26.79  on 2 df,   p=2e-06
coxFit<-coxph(Surv(coxDataRL$TimeMeasure, coxDataRL$Event == 1)~coxDataRL$rltl_change + coxDataRL$abs_rltl_change+ coxDataRL$average_resid+ as.numeric(paste(coxDataRL$amp_305_lactation)))

summary(coxFit)

See if relationship remains statistically significant when only animals with more than 2 samples are included in the analysis

coxData2<-data.frame(recoded_id=u_red_more_2$recoded_id, dob = u_red_more_2$dob, herd_life = u_red_more_2$herd_life, cull_date = u_red_more_2$cull_date, average_resid_M2 = u_red_more_2$average_resid_M2, rltl_change_M2 = u_red_more_2$rltl_change_M2, abs_rltl_change_M2 = u_red_more_2$abs_rltl_change_M2, amp_305_lactation = u_red_more_2$amp_305_lactation, av_rltl_1_res = u_red_more_2$av_rltl_1_res)




coxData2$CutoffDate<-"01/01/2017"
coxData2$Event<-ifelse(coxData2$cull_date=="NULL", 0, 1)
coxData2$TimeMeasure<-ifelse(coxData2$Event==0, as.Date(coxData2$CutoffDate, "%d/%m/%Y")- as.Date(coxData2$dob, "%d/%m/%Y"), as.numeric(paste(coxData2$herd_life)))
## Warning in ifelse(coxData2$Event == 0, as.Date(coxData2$CutoffDate, "%d/%m/%Y")
## - : NAs introduced by coercion
coxData2$average_resid_3G<-ntile(coxData2$average_resid_M2,3)
coxData2$rltl_change_3G<-ntile(coxData2$rltl_change_M2,3)
coxData2$abs_rltl_change_3G<-ntile(coxData2$abs_rltl_change_M2,3)

# models:
##Average lifetime telomere length
### condinuous variable:
coxFit<-coxph(Surv(coxData2$TimeMeasure, coxData2$Event == 1)~coxData2$average_resid_M2)

summary(coxFit)
## Call:
## coxph(formula = Surv(coxData2$TimeMeasure, coxData2$Event == 
##     1) ~ coxData2$average_resid_M2)
## 
##   n= 277, number of events= 213 
## 
##                              coef exp(coef) se(coef)      z Pr(>|z|)
## coxData2$average_resid_M2 -0.5990    0.5494   0.6041 -0.992    0.321
## 
##                           exp(coef) exp(-coef) lower .95 upper .95
## coxData2$average_resid_M2    0.5494       1.82    0.1682     1.795
## 
## Concordance= 0.524  (se = 0.022 )
## Likelihood ratio test= 0.99  on 1 df,   p=0.3
## Wald test            = 0.98  on 1 df,   p=0.3
## Score (logrank) test = 0.98  on 1 df,   p=0.3
### discrete scale

coxFit<-coxph(Surv(coxData2$TimeMeasure, coxData2$Event == 1)~coxData2$average_resid_3G)
summary(coxFit)
## Call:
## coxph(formula = Surv(coxData2$TimeMeasure, coxData2$Event == 
##     1) ~ coxData2$average_resid_3G)
## 
##   n= 277, number of events= 213 
## 
##                               coef exp(coef) se(coef)      z Pr(>|z|)
## coxData2$average_resid_3G -0.10350   0.90168  0.08528 -1.214    0.225
## 
##                           exp(coef) exp(-coef) lower .95 upper .95
## coxData2$average_resid_3G    0.9017      1.109    0.7629     1.066
## 
## Concordance= 0.526  (se = 0.021 )
## Likelihood ratio test= 1.47  on 1 df,   p=0.2
## Wald test            = 1.47  on 1 df,   p=0.2
## Score (logrank) test = 1.48  on 1 df,   p=0.2
coxFit<-survfit(Surv(coxData2$TimeMeasure, coxData2$Event == 1)~coxData2$average_resid_3G)
#### plot results for discrete scale

survp<- ggsurvplot(coxFit, data=coxData2, pval = TRUE, ggtheme = theme_minimal(), risk.table.y.text.col = T, risk.table.y.text = FALSE, legend.labs=c("Tertile 1 (shortest mean RLTL)", "Tertile 2", "Tertile 3 (longest mean RLTL)"))+theme_survminer(40)

survp

#ggsave(file = "/Users/luise/Documents/Work/CattleTelomereLength/AgingCell/Images/AverageRLTL.pdf", print(survp))

Mean RLTL change

coxFit<-coxph(Surv(coxData2$TimeMeasure, coxData2$Event == 1)~coxData2$rltl_change_M2)

summary(coxFit)
## Call:
## coxph(formula = Surv(coxData2$TimeMeasure, coxData2$Event == 
##     1) ~ coxData2$rltl_change_M2)
## 
##   n= 277, number of events= 213 
## 
##                             coef exp(coef) se(coef)      z Pr(>|z|)  
## coxData2$rltl_change_M2 -3.46163   0.03138  1.36042 -2.545   0.0109 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                         exp(coef) exp(-coef) lower .95 upper .95
## coxData2$rltl_change_M2   0.03138      31.87  0.002181    0.4515
## 
## Concordance= 0.544  (se = 0.023 )
## Likelihood ratio test= 6.42  on 1 df,   p=0.01
## Wald test            = 6.47  on 1 df,   p=0.01
## Score (logrank) test = 6.46  on 1 df,   p=0.01
### discrete scale

coxFit<-coxph(Surv(coxData2$TimeMeasure, coxData2$Event == 1)~coxData2$rltl_change_3G)
summary(coxFit)
## Call:
## coxph(formula = Surv(coxData2$TimeMeasure, coxData2$Event == 
##     1) ~ coxData2$rltl_change_3G)
## 
##   n= 277, number of events= 213 
## 
##                             coef exp(coef) se(coef)      z Pr(>|z|)
## coxData2$rltl_change_3G -0.14184   0.86776  0.08938 -1.587    0.113
## 
##                         exp(coef) exp(-coef) lower .95 upper .95
## coxData2$rltl_change_3G    0.8678      1.152    0.7283     1.034
## 
## Concordance= 0.534  (se = 0.021 )
## Likelihood ratio test= 2.52  on 1 df,   p=0.1
## Wald test            = 2.52  on 1 df,   p=0.1
## Score (logrank) test = 2.52  on 1 df,   p=0.1
coxFit<-survfit(Surv(coxData2$TimeMeasure, coxData2$Event == 1)~coxData2$rltl_change_3G)
#### plot results for discrete scale

survp<- ggsurvplot(coxFit, data=coxData2, pval = TRUE, ggtheme = theme_minimal(), risk.table.y.text.col = T, risk.table.y.text = FALSE, legend.labs=c("Tertile 1 (most RLTL shortening)", "Tertile 2", "Tertile 3 (least RLTL shortening)"))+theme_survminer(40)

survp

mean absolute RLTL change

coxFit<-coxph(Surv(coxData2$TimeMeasure, coxData2$Event == 1)~coxData2$abs_rltl_change_M2)

summary(coxFit)
## Call:
## coxph(formula = Surv(coxData2$TimeMeasure, coxData2$Event == 
##     1) ~ coxData2$abs_rltl_change_M2)
## 
##   n= 277, number of events= 213 
## 
##                              coef exp(coef) se(coef)     z Pr(>|z|)
## coxData2$abs_rltl_change_M2 1.522     4.579    1.159 1.312    0.189
## 
##                             exp(coef) exp(-coef) lower .95 upper .95
## coxData2$abs_rltl_change_M2     4.579     0.2184    0.4719     44.44
## 
## Concordance= 0.528  (se = 0.024 )
## Likelihood ratio test= 1.67  on 1 df,   p=0.2
## Wald test            = 1.72  on 1 df,   p=0.2
## Score (logrank) test = 1.72  on 1 df,   p=0.2
### discrete scale

coxFit<-coxph(Surv(coxData2$TimeMeasure, coxData2$Event == 1)~coxData2$abs_rltl_change_3G)
summary(coxFit)
## Call:
## coxph(formula = Surv(coxData2$TimeMeasure, coxData2$Event == 
##     1) ~ coxData2$abs_rltl_change_3G)
## 
##   n= 277, number of events= 213 
## 
##                                coef exp(coef) se(coef)     z Pr(>|z|)  
## coxData2$abs_rltl_change_3G 0.16969   1.18494  0.08624 1.968   0.0491 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                             exp(coef) exp(-coef) lower .95 upper .95
## coxData2$abs_rltl_change_3G     1.185     0.8439     1.001     1.403
## 
## Concordance= 0.537  (se = 0.021 )
## Likelihood ratio test= 3.88  on 1 df,   p=0.05
## Wald test            = 3.87  on 1 df,   p=0.05
## Score (logrank) test = 3.89  on 1 df,   p=0.05
coxFit<-survfit(Surv(coxData2$TimeMeasure, coxData2$Event == 1)~coxData2$abs_rltl_change_3G)
#### plot results for discrete scale

survp<- ggsurvplot(coxFit, data=coxData2, pval = TRUE, ggtheme = theme_minimal(), risk.table.y.text.col = T, risk.table.y.text = FALSE, legend.labs=c("Tertile 1 (least absolute RLTL change)", "Tertile 2", "Tertile 3 (most absolute RLTL change)"))+theme_survminer(40)

survp

Some more figures:

variable<-as.numeric(paste(data$sample_interval))
## Warning: NAs introduced by coercion
meanVar<-mean(variable, na.rm= TRUE)
sdVar<-sd(variable, na.rm=TRUE)

ggplot(data, aes(x = variable)) + 
    geom_histogram(colour= "black", fill = mycoloursP[7]) + xlab("Sample interval in days") +   geom_vline(xintercept= meanVar,  colour="red")+
  theme_classic(20)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning: Removed 305 rows containing non-finite values (stat_bin).

uData<-data[!duplicated(data$recoded_id),]

variable<-as.factor(paste(uData$birth_year))

ggplot(uData, aes(x = variable)) + 
    geom_bar(colour= "black", fill = mycoloursP[8]) + xlab("Birth year")+ theme_classic(20)

#Average milk production

variable<-as.numeric(paste(uData$amp_305_lactation))
## Warning: NAs introduced by coercion
meanVar<-mean(variable, na.rm= TRUE)
sdVar<-sd(variable, na.rm=TRUE)

ggplot(uData, aes(x = variable)) + 
    geom_histogram(colour= "black", fill = mycoloursP[9]) + xlab("Average milk production") +   geom_vline(xintercept= meanVar,  colour="red")+
  theme_classic(20)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning: Removed 52 rows containing non-finite values (stat_bin).

# cull reasons

variable<-as.factor(paste(uDeadData$cull_reason))

ggplot(uDeadData, aes(x = variable)) + 
    geom_bar(colour= "black", fill = mycoloursP[10]) + xlab("Cull reason") +
  theme_classic(10)+
  theme(axis.text.x = element_text(angle = 90, hjust = 1))

# fertility events
newdata <- subset(uData, uData$health_event_fertil_all != "NULL") 
variable<-as.integer(paste(newdata$health_event_fertil_all))

ggplot(newdata, aes(x = variable)) + 
    geom_bar(colour= "black", fill = mycoloursP[11]) + xlab("Fertility events")+
  theme_classic(20)+
  scale_x_continuous(breaks=c(0:12))

#mastitis events
newdata <- subset(uData, uData$health_event_mastitis_all != "NULL") 
variable<-as.integer(paste(newdata$health_event_mastitis_all))

ggplot(newdata, aes(x = variable)) + 
    geom_bar(colour= "black", fill = mycoloursP[12]) + xlab("Mastitis events")+
  theme_classic(20)+
  scale_x_continuous(breaks=c(0:8))

#lameness events
newdata <- subset(uData, uData$health_event_lame_all != "NULL") 
variable<-as.integer(paste(newdata$health_event_lame_all))

ggplot(newdata, aes(x = variable)) + 
    geom_bar(colour= "black", fill = mycoloursP[13]) + xlab("Lameness events")+
  theme_classic(20)+
  scale_x_continuous(breaks=c(0:11))

saving data

sessionInfo()
## R version 3.6.2 (2019-12-12)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS Catalina 10.15.5
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_GB.UTF-8/en_GB.UTF-8/en_GB.UTF-8/C/en_GB.UTF-8/en_GB.UTF-8
## 
## attached base packages:
## [1] grid      stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] ggsci_2.9         survminer_0.4.7   ggpubr_0.4.0      survival_3.2-3   
##  [5] dplyr_1.0.0       plyr_1.8.6        ggthemes_4.2.0    magrittr_1.5     
##  [9] stringr_1.4.0     data.table_1.12.8 mgcv_1.8-31       gridExtra_2.3    
## [13] lmerTest_3.1-2    lme4_1.1-23       Matrix_1.2-18     nlme_3.1-148     
## [17] ggplot2_3.3.2    
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_1.0.5          lattice_0.20-41     tidyr_1.1.0        
##  [4] zoo_1.8-8           digest_0.6.25       R6_2.4.1           
##  [7] cellranger_1.1.0    backports_1.1.8     evaluate_0.14      
## [10] highr_0.8           pillar_1.4.6        rlang_0.4.7        
## [13] curl_4.3            readxl_1.3.1        minqa_1.2.4        
## [16] car_3.0-8           nloptr_1.2.2.2      rmarkdown_2.3      
## [19] labeling_0.3        splines_3.6.2       statmod_1.4.34     
## [22] foreign_0.8-76      munsell_0.5.0       broom_0.7.0        
## [25] compiler_3.6.2      numDeriv_2016.8-1.1 xfun_0.15          
## [28] pkgconfig_2.0.3     htmltools_0.5.0     tidyselect_1.1.0   
## [31] tibble_3.0.3        km.ci_0.5-2         rio_0.5.16         
## [34] crayon_1.3.4        withr_2.2.0         MASS_7.3-51.6      
## [37] xtable_1.8-4        gtable_0.3.0        lifecycle_0.2.0    
## [40] KMsurv_0.1-5        scales_1.1.1        zip_2.0.4          
## [43] stringi_1.4.6       carData_3.0-4       farver_2.0.3       
## [46] ggsignif_0.6.0      ellipsis_0.3.1      survMisc_0.5.5     
## [49] generics_0.0.2      vctrs_0.3.1         boot_1.3-25        
## [52] openxlsx_4.1.5      tools_3.6.2         forcats_0.5.0      
## [55] glue_1.4.1          purrr_0.3.4         hms_0.5.3          
## [58] abind_1.4-5         yaml_2.2.1          colorspace_1.4-1   
## [61] rstatix_0.6.0       knitr_1.29          haven_2.3.1